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Abstract 

A  new  adaptive  split-domain  harmonic  balance  computational  fluid  dynamics 
(CFD)  method  is  developed  to  solve  highly  nonlinear  time-periodic  flows  such  as 
those  found  in  transonic  turbomachinery.  The  basic  harmonic  balance  CFD  method 
transforms  an  unsteady  time-periodic  problem  into  a  steady-state  problem  by  as¬ 
suming  a  solution  in  the  form  of  a  Fourier  series  in  time.  The  new  method  employs 
a  unique  multi-domain  split-operator  solution  technique  to  remove  a  large-series  sta¬ 
bility  restriction  present  in  previous  harmonic  balance  CFD  approaches.  The  new 
method  also  minimizes  the  computational  work  required  to  obtain  a  harmonic  bal¬ 
ance  solution  by  adapting  the  frequency  content  to  the  flow,  starting  with  a  small 
number  of  Fourier  frequencies  and  augmenting  the  frequency  content  in  each  cell  as 
necessary  to  capture  local  flow  physics.  The  method  reduces  compute  times  by  allow¬ 
ing  larger  integration  time  steps,  reducing  Fourier  transform  calculations,  and  reduc¬ 
ing  overall  problem  size.  Split-domain  solutions  to  the  1-D  inviscid  Burgers’  equation 
are  computed  with  up  to  97  frequencies,  demonstrating  improved  stability.  Differ¬ 
ences  between  harmonic  balance  solutions  and  time-accurate  solutions  are  found  to 
be  asymptotic  with  Fourier  series  length.  The  adaptive  split-domain  approach  is 
applied  to  the  1-D  and  quasi-l-D  Euler  equations.  Supersonic  and  subsonic  Euler 
calculations  show  that  the  adapted  and  non-adapted  harmonic  balance  solutions  are 
equivalent.  Accurate  adapted  quasi-l-D  Euler  solutions  for  a  supersonic/subsonic 
diverging  nozzle  with  periodic  unsteady  outflow  conditions  are  generated  in  86% 
less  time  than  an  equivalent  non-adapted  split-domain  solution,  demonstrating  the 
benefit  of  adapting  frequency  content  to  local  flow  conditions. 
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ADAPTIVE  HARMONIC  BALANCE  METHOD 
EOR  UNSTEADY,  NONLINEAR, 
ONE-DIMENSIONAL  PERIODIC  ELOWS 


/.  Introduction 

High-fidelity  numerical  simulations  of  fluid  flow  through  transonic  turboma¬ 
chinery  are  of  considerable  interest  to  designers  of  compressors  and  turbines  in  mod¬ 
ern  jet  engines.  Solutions  can  be  obtained  with  conventional  time-accurate  compu¬ 
tational  fluid  dynamics  (CFD)  codes,  but  the  considerable  time  required  to  generate 
these  solutions  limits  their  utility  to  the  designer.  For  the  class  of  problems  where 
the  flow  can  be  assumed  to  be  fully  developed  and  periodic  in  time,  such  as  flow 
past  a  rotor  with  oscillating  blades  or  flow  through  a  rotor-stator,  time-accurate 
calculations  can  be  particularly  inefficient.  It  is  usually  necessary  to  step  through 
many  disturbance  periods  before  a  fully  developed  solution  is  reached.  To  achieve 
shorter  computation  times  for  this  class  of  problem,  CFD  techniques  have  been  de¬ 
veloped  that  take  advantage  of  the  time-periodic  nature  of  the  flow.  These  include 
the  time-linearization  technique  (1,  2,  3,  4,  5,  6,  7,  8,  9,  10,  11,  12),  the  time-averaging 
technique  (13,  14,  15),  and  the  harmonic  balance  technique  (16,  17,  18,  19). 

Of  the  three  CFD  techniques  developed  specihcally  for  time-periodic  flows,  the 
harmonic  balance  technique  is  most  suitable  for  modeling  the  large  disturbances  and 
strong  nonlinearities  found  in  transonic  turbomachinery.  All  three  techniques  are 
closely  related  in  that  they  all  assume  a  harmonic  form  for  the  unsteady  flow,  recast 
the  unsteady  problem  as  a  steady-state  calculation,  and  employ  convergence  accel¬ 
eration  techniques  to  reduce  computation  time.  However,  the  time-linearization  and 
time-averaging  techniques  have  restrictions  on  the  unsteadiness  that  are  not  present 
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in  the  harmonic  balance  method.  The  time-linearized  technique  models  unsteadiness 
as  a  small-amplitude  perturbation  linearized  about  an  nonlinear  steady-state  back¬ 
ground  flow.  Linearization  of  the  unsteadiness,  combined  with  a  small-disturbance 
assumption,  limits  applicability  to  transonic  flows.  In  the  time-averaged  technique 
a  small  harmonic  perturbation  is  time-averaged  with  the  steady-state  background 
equations  in  a  process  similar  to  Reynolds  averaging.  Additional  stress  terms  re¬ 
sult  that  capture  some  nonlinearities  in  the  unsteady  flow,  but  a  small  disturbance 
assumption  is  still  made.  Finally,  both  the  time-linearized  and  time-averaged  tech¬ 
niques  solve  for  only  one  harmonic  at  a  time.  Multiple-harmonic  solutions  are  built 
up  by  superposing  several  single-harmonic  solutions,  losing  the  effect  of  nonlinear 
harmonic  coupling.  The  harmonic  balance  method,  on  the  other  hand,  has  no  small 
disturbance  restriction,  and  solves  a  set  of  nonlinear  equations  for  several  harmonics 
simultaneously,  more  accurately  capturing  coupled  nonlinear  effects.  This  makes  it 
uniquely  suited  for  calculating  flow  through  transonic  turbomachinery. 

The  fidelity  of  a  harmonic  balance  solution  is  dependent  on  grid  density,  on 
the  number  of  harmonics  included,  and  on  the  flow  being  modeled.  On  a  given  grid, 
a  flow  that  is  smoothly  unsteady  (i.e.,  without  moving  discontinuities)  will  require 
fewer  harmonics  than  a  flow  containing  a  moving  shock  to  achieve  the  same  level  of 
fidelity.  Because  the  computational  cost  increases  with  each  harmonic  included  in 
the  solution,  it  is  desirable  to  use  the  minimum  number  of  harmonics  needed  for  a 
given  problem. 

In  a  typical  transonic  turbomachinery  problem,  the  nature  of  the  flow  can 
vary  significantly  throughout  the  domain  of  interest.  This  is  illustrated  in  Fig.  1.1 
(20),  which  shows  experimental  time-pressure  plots  at  two  locations  on  a  single  inlet 
guide  vane  upstream  of  a  transonic  rotor.  Existing  harmonic  balance  implementa¬ 
tions  solve  for  a  constant  number  of  harmonics  over  the  entire  computational  domain, 
so  problems  containing  both  smooth  and  discontinuous  unsteadiness  require  a  com¬ 
promise  between  reduced  run  time  (fewer  harmonics),  and  fidelity  in  the  regions  of 
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95%  ClKid 


reproduced  by  permission  of  authors  (20) 


Figure  1.1  Experimental  Pressure  Data  on  Inlet  Guide  Vane  Upstream  of  Tran¬ 
sonic  Rotating  Compressor  Blade  Row  for  5%  Span  (top),  and  95% 
Span  (bottom) 

strong  nonlinearity.  The  goal  of  this  research  is  to  remove  the  need  for  compromise 
by  extending  the  harmonic  balance  technique  to  allow  a  variable  number  of  included 
harmonics,  and  to  automatically  identify  and  apply,  on  a  cell-by-cell  basis,  the  min¬ 
imum  number  of  harmonics  required  to  achieve  a  consistent  hdelity  throughout  the 
computational  domain. 

1.1  Overview  of  Harmonie  Balanee  Method 

The  harmonic  balance  method  has  been  used  for  many  years  as  a  means  of  an¬ 
alyzing  the  behavior  of  harmonic  ordinary  differential  equations  (ODEs)  (21).  The 
technique  consists  of  assuming  a  solution  in  the  form  of  a  truncated  Fourier  series 
with  a  predetermined  number  of  harmonics,  substituting  the  assumed  solution  into 
the  ODE,  and  algebraically  manipulating  the  results  to  collect  terms  of  like  fre¬ 
quency.  Any  resulting  terms  with  a  frequency  not  in  the  Fourier  series  are  dropped. 
Each  harmonic  is  then  balanced  by  requiring  that  like-frequency  terms  on  each  side 
of  the  equation  satisfy  equality  independently.  Balancing  results  in  a  system  of  cou- 


1-3 


pled  algebraic  equations  which  are  solved  for  the  Fourier  coefficients  of  the  assumed 
solution. 


When  the  same  technique  is  applied  to  a  partial  differential  equation  (PDF), 
the  result  is  a  system  of  ODEs  or  PDEs  that  are  solved  for  the  Fourier  coefficients. 
For  example,  consider  a  generic  scalar  one-dimensional  Conservation  Law  equation 
in  differential  form 


dt  dx 


(1.1) 


In  this  equation,  ^{x,t)  is  the  scalar  dependent  variable,  and  <F  is  a  flux  term  that 
depends  on 


If  the  time  response  of  the  conserved  variable  ^  at  an  arbitrary  point  in  space, 
X,  is  assumed  to  be  periodic  in  time  with  radian  frequency  u,  then  that  response  can 
be  approximated  by  a  truncated  complex  Fourier  series 


N 

ax,t)^  Cnix)e^^^\  (1.2) 

n=-N 


In  this  series,  the  coefficients,  c„,  are  functions  of  x  only,  while  the  exponential  terms 
are  functions  of  time  only.  To  derive  the  harmonic  balance  form  of  the  conserva¬ 
tion  equation,  the  approximating  series  is  substituted  into  the  conservation  equation 
which,  after  some  algebraic  manipulation,  is  transformed  into  a  system  of  ordinary 
differential  equations  dependent  only  on  x,  which  are  solved  for  the  unknown  Fourier 
coefficients.  Because  the  variable  ^{x,t)  is  real,  the  Fourier  coefficients  correspond¬ 
ing  to  ±n  are  complex  conjugates  of  each  other.  Equations  associated  with  negative 
frequencies  can  therefore  be  dropped,  leaving  iV  -|-  1  complex  differential  equations 
for  2N  +  1  unknowns.  In  vector  form,  these  equations  are  given  by 

^  +  S(0  =  0  (1.3) 
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where  ^  is  the  vector  of  Fourier  coefficients,  c„,  <F  is  a  vector  of  harmonic  balance  flux 
terms,  and  S  is  a  vector  of  source  term  arising  from  the  time  derivative  in  Eq.  1.1. 

For  problems  of  interest,  Eq.  1.3  is  highly  nonlinear.  A  pseudo-time  derivative 
of  the  dependent  vector  ^  is  added  (Eq.  1.4)  so  that  a  time- marching  solution  method 
can  be  applied  to  hnd  the  steady-state  solution 


dr  dx 


+  S(0  =  0. 


(1.4) 


When  steady  state  is  reached,  the  added  pseudo-time  derivative  is  equal  to  zero,  and 
Eq.  1.3  is  recovered. 

Once  a  steady-state  solution  is  obtained,  useful  information  is  recovered  from 
the  calculated  Fourier  coefficients.  The  values  of  the  dependent  variable  at  any  point 
in  time  are  easily  reconstructed  from  the  computed  coefficients  by  performing  the 
summation  in  Eq.  1.2.  The  time-average  values  of  the  dependent  variable  are  also 
readily  obtained,  as  they  are  the  computed  coefficients  cq. 


Previous  Work  in  Harmonic  Balance.  The  hrst  use  of  the  harmonic  bal¬ 
ance  method  in  CFD  was  by  Hall,  Thomas,  and  Clark  (16),  who  implemented  a 
harmonic  balance  Reynolds-Averaged  Navier-Stokes  solver.  They  applied  it  to  a 
single  two-dimensional  compressor  blade  row  undergoing  forced  periodic  vibration 
under  transonic  flow  conditions.  This  configuration  contained  moderately  strong 
shocks,  but  limited  shock  motion.  Solutions  containing  up  to  7  frequencies  in  the 
approximating  series  were  calculated.  Solutions  containing  3-5  frequencies  showed 
good  agreement  with  standard  time-accurate  calculations  for  the  time-average  and 
hrst  harmonic  terms,  but  took  approximately  one-tenth  as  long  to  compute.  Solu¬ 
tions  with  7  frequencies  failed  to  converge.  In  a  follow-on  effort,  Thomas,  Dowell, 
and  Hall  (18)  coupled  an  inviscid  harmonic  balance  solver  with  a  linear  structural 
model  to  study  limit  cycle  oscillations  of  an  aeroelastic  system. 
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Additional  work  in  harmonic  balance  CFD  was  conducted  by  McMullen,  Jame¬ 
son,  and  Alonso  (17,  19),  who  investigated  several  aspects  of  the  harmonic  balance 
method  not  previously  addressed.  Using  a  somewhat  different  formulation  than  Hall 
et  ai,  they  performed  a  stability  analysis  of  the  method,  concluding  that  stability 
could  be  an  issue  when  large  numbers  of  frequencies  are  included  in  the  solution. 
They  applied  the  method  to  several  test  problems,  including  unsteady  1-D  chan¬ 
nel  flow,  2-D  oscillating  flow  behind  a  cylinder  in  crossflow,  and  a  pitching  airfoil. 
Their  results  show  good  agreement  with  analytical  and  numerical  predictions  for 
both  time  dependent  and  time-averaged  data.  Early  cylinder  crossflow  calculations 
hxed  the  fundamental  frequency  at  the  theoretical  frequency  of  vortex  shedding  (17). 
Later  cylinder  crossflow  and  pitching  airfoil  calculations  employed  a  gradient-based 
optimization  approach  to  dynamically  determine  the  correct  fundamental  frequency 
(19). 

1.2  Scope 

The  objective  of  this  research  is  to  extend  the  harmonic  balance  technique 
to  accurately  resolve  unsteady,  nonlinear,  periodic  flows  while  minimizing  compu¬ 
tational  cost.  To  achieve  this  objective,  the  solution  technique  must  be  able  to 
robustly  solve  for  as  many  harmonics  as  necessary  to  attain  the  desired  degree  of 
solution  hdelity.  To  maintain  overall  efficiency  and  minimize  computational  cost, 
higher  harmonics  should  only  be  included  where  required  by  local  flow  conditions. 
The  scope  of  this  research  is  dehned  by  the  following  thesis  statement  and  assump¬ 
tions: 


1.2.1  Thesis  Statement.  A  spatially-adaptive  harmonic  balance  method 
can  he  implemented  to  accurately  and  efficiently  compute  a  stationary  time-periodic 
flow  field  containing  regions  of  smooth  flow  and  regions  with  strong  moving  discon¬ 
tinuities  by  automatically  adjusting  the  number  of  freguencies  in  the  solution,  on  a 
point-by-point  basis,  to  match  local  flow  conditions. 
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1.2.2  Assumptions.  The  implementation  and  effectiveness  of  the  spatially- 
adaptive  harmonic  balance  technique  are  not  strongly  dependent  on  the  spatial  di¬ 
mension  of  the  problem,  and  will  be  adequately  demonstrated  in  one  spatial  dimen¬ 
sion.  Efficiency  will  be  measured  by  comparing  average  frequency  content  and  run 
time  with  a  comparable  non-adapted  harmonic  balance  implementation.  Accuracy 
will  be  determined  through  comparison  with  conventional  time-accurate  results  and 
non-adapted  harmonic  balance  results. 

1.3  Research  Approach 

Several  one-dimensional  CFD  codes  were  developed  to  investigate  different 
aspects  of  the  adaptive  harmonic  balance  method.  The  initial  work  investigated 
harmonic  balance  solutions  of  the  inviscid  Burgers’  equation,  a  scalar  analog  of  the 
Euler  equation,  subject  to  a  variety  of  periodic  boundary  conditions.  Next,  subsonic 
and  supersonic  harmonic  balance  solutions  to  Euler’s  equation  were  investigated. 
Finally,  the  adaptive  split-domain  harmonic  balance  method  was  demonstrated  by 
solving  the  quasi-l-D  Euler’s  equation  for  an  unsteady  diverging  nozzle  conhguration. 

Inviscid  Burgers’  Equation:  Different  implementations  of  the  harmonic  bal¬ 
ance  Burgers’  equation  were  tested  to  determine  which  implementation  provided  the 
best  performance,  robustness,  and  accuracy.  A  new  split-domain  implementation 
was  shown  to  be  a  superior  harmonic  balance  implementation  when  large  numbers 
of  harmonics  were  included.  The  impact  of  shock  strength,  fundamental  frequency, 
and  grid  density  on  the  harmonic  balance  solution  were  investigated. 

Euler’s  Equation:  An  adaptive  split-domain  harmonic  balance  Euler  solver  was 
written  to  investigate  the  frequency  augmentation  approach  and  to  develop  methods 
for  properly  treating  transitions  from  one  frequency  to  another  at  arbitrary  points 
in  the  computational  grid.  A  study  of  the  impact  of  user-specified  parameters  on 
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the  performance  of  the  adaptive  algorithm  for  both  supersonic  and  subsonic  flows 
was  conducted. 


Quasi- 1-D  Euler’s  Equation:  The  adaptive  split-domain  harmonic  balance  al¬ 
gorithm  is  demonstrated  by  solving  unsteady  diverging  nozzle  flow.  This  conhgu- 
ration  is  representative  of  anticipated  production  applications,  containing  regions  of 
both  low  and  high  frequency  content,  with  abrupt  transitions  between  the  two. 


1.4  Document  Organization 

The  remainder  of  this  document  is  organized  as  follows: 

Chapter  H:  Details  the  theory  and  implementation  of  the  adaptive  split-domain 
harmonic  balance  method. 

Chapter  IH:  Records  the  analysis  and  results  of  the  harmonic  balance  method 
applied  to  the  inviscid  1-D  Burgers  equation  for  a  variety  of  periodic  boundary 
conditions. 

Chapter  IV:  Records  the  analysis  and  results  of  the  adaptive  split-domain 
harmonic  balance  method  applied  to  the  one-dimensional  Euler  equation  for  both 
supersonic  and  subsonic  boundary  conditions. 

Chapter  V:  Presents  a  demonstration  of  the  adaptive  split-domain  harmonic 
balance  technique  applied  to  an  unsteady  quasi- 1-D  nozzle  flow. 

Chapter  VI:  Summarizes  the  results  and  conclusions  of  the  current  research 

Chapter  VII:  Discusses  future  research  topics  suggested  by  the  current  re¬ 
search. 
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11.  Theory  and  Implementation  of  Adaptive  Split-Domain  Harmonic 

Balance 

2. 1  Introduction 

In  this  chapter,  a  theoretical  and  practical  foundation  for  the  adaptive  split- 
domain  harmonic  balance  CFD  method  is  established.  The  chapter  begins  by  ex¬ 
panding  on  the  brief  overview  of  harmonic  balance  in  CFD  presented  in  Section  1.1 
with  a  more  in-depth  derivation  of  the  basic  harmonic  balance  method  and  its  appli¬ 
cation  to  the  solution  of  conservation  equations.  This  is  followed  by  a  brief  discussion 
of  previous  harmonic  balance  CFD  implementations,  and  a  detailed  description  of  the 
new  split-domain  harmonic  balance  approach.  The  focus  of  the  discussion  then  turns 
to  frequency  augmentation  and  the  adaptive  harmonic  balance  algorithm.  Finally, 
the  chapter  concludes  with  a  discussion  of  the  implementation  of  full  approxima¬ 
tion  scheme  (FAS)  multigrid  convergence  acceleration  with  adaptive  split-domain 
harmonic  balance. 

The  discussion  that  follows  is  general.  No  specihc  solver  implementation  is 
described,  and  no  assumptions  are  made  about  the  specific  conservation  equations 
being  solved  or  the  discretization  schemes  used  to  obtain  a  numerical  solution.  Be¬ 
cause  of  this,  the  concepts  and  algorithms  described  below  are  easily  adapted  to  new 
applications.  In  the  chapters  that  follow,  some  or  all  of  the  ideas  discussed  below  are 
applied  to  a  MacCormack  discretization  of  the  1-D  Burgers’  equation  (Chapter  III), 
as  well  as  cell-centered  finite  volume  formulations  of  the  1-D  Euler’s  equation  (Chap¬ 
ter  IV)  and  quasi-l-D  Euler’s  equation  (Chapter  V). 

2.2  Derivation  of  the  Harmonic  Balance  Equations 

In  Section  1.1,  a  brief  overview  of  the  application  of  the  harmonic  balance 
method  to  a  scalar  conservation  law  was  presented.  The  following  discussion  expands 
on  that  brief  overview,  providing  a  detailed  derivation  of  the  harmonic  balance  form 
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of  the  conservation  law,  and  covering  some  aspects  of  numerical  implementation. 
The  derivation  is  carried  out  based  on  the  generic  scalar  conservation  law,  Eq.  1.1, 
repeated  here  as  Eq.  2.1 


,  dm  n 

dt  dx 


(2.1) 


A  similar  derivation  for  a  specihc  conservation  law,  the  1-D  inviscid  Burgers’  equation 
(Eq.  3.1)  is  presented  in  Appendix  A. 


As  discussed  in  Section  1.1,  the  harmonic  balance  method  centers  on  the  as¬ 
sumption  that  the  dependent  variable,  can  be  approximated  by  a  truncated  Fourier 
series 

N 

Y.  (2.2) 

n=-N 

In  this  series,  N  is  the  number  of  positive  and  negative  frequencies  in  the  truncated 
Fourier  series,  and  u  is  the  fundamental  radian  frequency  of  the  series,  equal  to  27r/P 
where  P  is  the  period  of  oscillation  of  the  physical  system.  The  complex  coefficients, 
Cn,  are  functions  of  the  spatial  variable  only,  while  the  exponentials  are  functions  of 
time.  Substituting  Eq.  2.2  into  the  hrst  term  on  the  left  hand  side  of  Eq.  2.1  and 
evaluating  the  time  derivative  yields 


dt 


N 

n=—N 


(2.3) 


This  expression  is  the  origin  of  the  harmonic  balance  source  term. 

A  similar  substitution  is  performed  for  the  second  term  on  the  left  hand  side 
of  Eq.  2.1.  In  this  case,  substitution  requires  the  evaluation  of  the  flux  function,  <h, 
with  the  approximating  series  as  an  argument.  Upon  evaluation,  the  resulting  flux 
is  algebraically  simplihed  and  terms  of  like  frequency  are  collected,  resulting  in  an 
expression  of  the  form 

M 

~  (2.4) 

n=—M 
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where  (pn  contains  the  sum  of  all  coefficients  of  the  exponential  If  the  evaluation 
of  <h  requires  products  or  integer  powers  of  the  dependent  variable  then  M  >  N . 
As  part  of  the  harmonic  balance  approximation,  all  exponential  terms  with  frequency 
greater  than  iiVci;  are  discarded,  leaving  only  terms  with  frequencies  present  in  the 
approximating  series 

N 

'!>(«)  »  Y.  P.5) 

n=-N 

As  a  simple  example  of  the  derivation  of  a  harmonic  balance  flux,  consider  a 
flux  function  $(.0  =  and  an  approximating  Fourier  series  with  iV  =  1.  Evaluation 
of  the  flux  with  the  approximating  series  as  the  argument  yields 

.  (2.6) 

Expanding  the  right-hand  side  of  Eq.  2.6  gives 

c_]^e  -|“  c_iCoe  -h  c_iCi  cqC—ic — iujt  (2.7) 

cl  +  cocie'"*  +  cic_i  +  cicoe'"*  + 

+  2c_iCoe  +  2c_iCi  +  Cq  + 

2coCie*‘^‘ +  (2.8) 

Equation  2.8  contains  exponential  terms  with  frequencies  ±2cc;.  Since  the  assumed 
solution  contains  frequencies  only  up  to  icu,  these  higher  frequency  terms  are  dis¬ 
carded.  The  remaining  terms  in  Eq.  2.8  are  the  harmonic  balance  flux  terms  for  this 
simple  example: 


0-1  —  2c_iCo 


00  —  2c_lCl  -|-  Cq 


0_1  =  2ciCo.  (2.9) 
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Substitution  of  Eq.  2.3  and  Eq.  2.5  into  the  left  hand  side  of  Eq.  2.1  and 
collection  of  terms  with  like  exponentials  yields  the  series  approximation  form  of  the 
conservation  equation 

N 

n=-N 

At  this  point,  the  harmonic  terms  are  required  to  balance  across  the  equality;  each 
frequency  is  required  to  satisfy  equality  independently  of  the  other  frequencies.  Since 
the  right  hand  side  of  Eq.  2.10  is  identically  zero,  each  of  the  terms  in  the  summation 
on  the  left  hand  side  must  be  equal  to  zero,  i.e.. 


inujCr,  + 


cl(f)n 

dx 


(2.10) 


tnojc^i  ~\~  — ; — 
dx 


^incot 


0 


-N  <n<N.  (2.11) 


The  terms  inside  the  square  brackets  must  be  identically  equal  to  zero  for  each  of  the 
2N  +  1  complex  equations  represented  by  Eq.  2.11  to  hold  for  all  time.  The  expo¬ 
nential  terms  are  dropped,  leaving  a  system  of  2N  +  1  coupled  ordinary  differential 
equations  for  2N  +  1  complex  coefficients.  This  system  can  be  further  simplihed 
because  the  dependent  variable,  ^{x,t),  is  real  and  periodic.  The  negative  frequency 
Fourier  coefficients  are  thus  the  complex  conjugates  of  the  corresponding  positive 
frequency  coefficients,  i.e.,  C-n  =  Cn-  Furthermore,  since  the  zero-frequency  term  cq 
must  be  real,  the  total  number  of  unknowns  is  2N  +  1:  the  real  and  imaginary  parts 
of  N  positive-frequency  coefficients,  plus  the  real  zero-frequency  coefficient.  Since 
the  coefficients  of  the  negative  frequencies  are  not  independent,  the  corresponding 
equations  need  not  be  solved.  Keeping  just  the  positive  frequency  equations  and 
eliminating  the  exponential  terms,  Eq.  2.11  reduces  to 


d(j)n 

tnUlCn  H - ; — 

dx 


0 


0  <  n  <  iV,  (2.12) 


or  in  vector  form. 


dm 

dx 


+m  =  o 


(2.13) 
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with 


Co 

0 

00 

Cl 

S(0  = 

iuci 

m  = 

01 

Cn 

iNojcjsi 

4>n 

(2,14) 


To  solve  Eq.  2.13,  a  pseudo-time  derivative  of  the  vector  of  Fourier  coefficients,  is 
added 


(2.15) 


A  time-marching  scheme  is  used  to  drive  the  solution  to  steady  state.  Once  steady 
state  is  reached,  the  pseudo-time  derivative  vanishes,  and  Eq.  2.13  is  recovered. 


Eq.  2.15  is  the  vector  harmonic  balance  form  of  the  scalar  one-dimensional 
conservation  equation  (Eq.  2.1).  Derivation  of  the  harmonic  balance  form  for  multi¬ 
dimensional  and  vector  conservation  equations  is  similar  to  that  of  the  scalar  equa¬ 
tion.  For  a  multi-dimensional  problem,  multiple  flux  derivatives  must  be  computed, 
but  the  basic  harmonic  balance  form  remains  unchanged.  In  the  case  of  a  vector 
equation  (e.g.  the  Euler  equation),  each  of  the  conserved  (or  primitive)  variables 
is  expanded  in  a  separate  Fourier  series.  The  harmonic  balance  solution  vector  is 
a  concatenation  of  the  Fourier  coefficients  for  all  the  approximated  variables.  An 
example  of  a  vector  equation  derivation  is  found  in  (16). 

Evaluating  the  harmonic  balance  flux  vector,  <l>,  can  be  computationally  expen¬ 
sive.  For  a  simple  flux  like  that  of  the  scalar  1-D  inviscid  Burgers’  equation  (Eq.  3.1), 
the  asymptotic  complexity  of  the  harmonic  balance  flux  calculation  is  (9(iV^).  For 
a  more  complex  vector  equation  such  as  the  Euler  equation  (Eq.  4.1),  the  asymp¬ 
totic  complexity  lies  between  0{N^)  and  O(iV^)  (16).  Considerable  computational 
savings  can  be  realized  by  taking  a  multi-domain  approach  for  the  flux  calculation. 
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A  multi-domain  approach  takes  advantage  of  the  relationship  between  the  com¬ 
puted  Fourier  coefficients,  and  their  inverse  discrete  Fourier  transform  (IDFT) 
(16).  Given  the  vector  of  Fourier  coefficients  and  JF,  a  discrete  Fourier  transform 
(DFT)  operator  that  produces  positive-frequency  Fourier  coefficients  from  a  set  of 
real  numbers,  dehne  ^  such  that 

e  =  H  (2.16a) 

e  =  (2.16b) 

The  vector  is  a  real  vector  of  length  2N  +  1  that  contains  values  of  the  dependent 
variable,  sampled  at  times  t  =  (0,  At,  2Af,  •  •  ■  ,2NAt),  with  At  =  ^(2^+1)  ’ 

the  period  of  oscillation  divided  by  the  number  of  samples.  A  similar  relationship  is 
assumed  for  the  flux  vector,  $,  and  its  IDFT  (16). 

Given  these  relationships,  the  harmonic  balance  flux,  <h,  can  be  closely  approx¬ 
imated  by  calculating  ^  from  applying  the  time  domain  flux  function  <F  to  each 
element  of  ^  to  obtain  a  vector  of  time-sampled  fluxes,  <h,  and  transforming  the  result 
back  to  the  frequency  domain  (16).  The  asymptotic  complexity  of  the  DFT/IDFT 
for  arbitrary  numbers  of  terms  is  O(iV^)  (22).  The  complexity  of  the  time-domain 
flux  calculation  is  0{N),  so  the  complexity  of  the  combined  operation  is  O(iV^). 
To  put  the  potential  computational  savings  in  perspective,  consider  an  Euler  flux 
calculation  with  N  =  45,  and  assume  a  computational  cost  proportional  to  for 
the  multi-domain  approach,  and  for  the  direct  approach.  Disregarding  any  con¬ 
stants  of  proportionality,  the  direct  approach  is  over  300  times  more  expensive  than 
the  multi-domain  approach. 

In  practice,  the  theoretical  computational  cost  savings  could  be  much  higher. 
The  asymptotic  complexity  of  a  modern  DFT/IDFT  algorithm  has  a  lower  bound 
of  D(Alog2  N)  (22).  The  lower  bound  is  typically  obtainable  only  when  the  number 
of  terms  being  transformed  is  a  power  of  2,  however.  Since  the  number  of  terms 
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transformed  in  the  harmonic  balance  calculation  is  always  odd,  the  ideal  performance 
is  never  achieved.  However,  very  good  performance  can  still  be  achieved  if  the 
number  of  terms  can  be  factored  into  integer  powers  of  small  primes  (22),  i.e.,  when 
2iV  + 1  =  2“3^5'^7‘^1H  for  some  integers  a,  b,  c,  d,  and  e.  So  the  actual  computational 
cost  will  lie  between  iVlog2  N  and  iV^,  depending  on  the  value  of  N.  For  the  above 
example  with  iV  =  45,  a  direct  flux  calculation  is  between  300  and  2400  times  more 
costly  than  a  multi-domain  calculation,  disregarding  constants  of  proportionality. 


2.3  Prior  Implementations  of  Harmonic  Balance  for  CFD 

Existing  harmonic  balance  CFD  solvers  have  taken  two  different  approaches  to 
implementing  the  harmonic  balance  equations.  Both  are  multi-domain  approaches. 
The  hrst  approach,  developed  by  Hall,  Thomas,  and  Clark  (16),  results  in  a  system 
of  equations  that  is  integrated  in  the  time  domain.  The  second,  developed  by  Mc¬ 
Mullen,  Jameson,  and  Alonso  (17),  results  in  a  system  of  equations  integrated  in  the 
frequency  domain. 

Derivation  of  the  time-domain  approach  begins  by  substituting  Eq.  2.16a  and 
the  multi-domain  flux  calculation  into  Eq.  2.15  to  obtain 


dr 


dx 


+  =  0. 


(2,17) 


An  IDFT  applied  to  Eq.  2.17  yields  the  hnal  form  of  the  time-domain  harmonic 
balance  approach. 

^  +  5®^  +  ^“^8(5^0  =  0.  (2.18) 

OT  OX 

Eq.  2.18  has  the  form  of  the  original  conservation  equation,  Eq.  2.1,  with  an  added 
source  term.  The  primary  advantage  of  this  formulation  is  that  an  existing  steady- 
state  solver  can  be  easily  modified  to  solve  the  harmonic  balance  equations. 
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Derivation  of  the  frequency-domain  formulation  begins  by  rewriting  Eq.  2.15 
as 

^  +  m  +  S(0  =  0,  (2.19) 

where  R{^)  is  a  residual  containing  the  flux  derivative  term,  (or  in  multiple 
dimensions,  the  sum  of  the  flux  derivatives  for  each  dimension)  as  well  as  any  terms 
required  by  a  specihc  discretization,  such  as  artihcial  dissipation  terms.  A  frequency 
domain/time  domain  relationship  is  assumed  for  the  residual,  such  that 

R{i)  =  (2.20) 

where  R{^)  contains  the  flux  derivatives  and  additional  terms  evaluated  at  2N  +  1 
points  in  time.  Substitution  of  Eq.  2.20  into  Eq.  2.19  yields  the  second  harmonic 
balance  form, 

^  +  .F:R(.F-iO  +  S(eH0.  (2.21) 

Despite  being  integrated  in  the  frequency  domain,  this  formulation  is  very  similar 
to  the  time-domain  formulation.  Due  to  the  linearity  of  T  and  if  the  same 

discretization  scheme  is  applied  to  both  formulations,  and  the  spatial  discretiza¬ 
tion  of  i?(0  contains  only  linear  difference  operators,  then  Eqs.  2.18  and  2.21  are 
mathematically  equivalent. 

2.J^  Split-Domain  Harmonic  Balance 

The  time-domain  and  frequency- domain  formulations  of  the  harmonic  balance 
equations  are  adequate  for  calculating  solutions  with  small  N ,  but  both  have  aspects 
that  reduce  their  utility  for  solutions  requiring  large  N.  The  primary  concern  with 
the  existing  formulations  is  the  stability  of  the  resulting  harmonic  balance  equations. 
A  stability  analysis  of  the  frequency-domain  approach  indicates  that  as  N  becomes 
large,  the  stability  limit  of  the  approach  becomes  restricted  (17),  requiring  a  smaller 
time  step  and  reducing  efficiency.  This  property  was  confirmed  experimentally  for 
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both  the  time-domain  and  frequency  domain  formulations,  as  well  as  for  the  direct 
formulation  (see  Section  3.4.2). 

The  second,  lesser  concern  with  the  existing  formulations  is  the  number  of 
transforms  required  per  point,  per  iteration.  A  multi-step  time  integration  technique 
such  as  an  explicit  3-,  4-,  or  5-stage  Runge-Kutta  scheme  requires  one  DFT/IDFT 
evaluation  per  point  for  each  stage  of  the  integrator.  For  small  N,  this  is  of  little 
consequence,  but  for  large  N  the  cost  of  multiple  DFT/IDFT  evaluations  could 
become  signihcant. 

A  new  split-domain  harmonic  balance  formulation  was  developed  that  ad¬ 
dresses  both  the  stability  concern  and  the  DFT/IDFT  evaluation  concern.  In  this 
formulation,  the  inhomogeneous  harmonic  balance  equation  (Eq.  2.15)  is  split  into  a 
homogeneous  partial  differential  equation  and  an  ordinary  differential  equation  (23). 


dji  dHii)  _ 

dr  dx 

^  +  =  0. 

dr 


(2.22a) 

(2.22b) 


An  approximate  solution  to  Eq.  2.15  is  obtained  by  sequentially  solving  Eqs.  2.22a 
and  2.22b,  using  the  solution  from  one  as  the  initial  condition  for  the  other.  For 
example,  one  could  solve  Eq.  2.22a  for  ^i{x,t)  subject  to  the  initial  condition 
.^i(a;,0)  =  .^(a;,0),  the  initial  condition  for  Eq.  2.15.  The  approximate  solution  is 
then  obtained  by  solving  Eq.  2.22b  for  ^2  with  the  initial  condition  .^2(0)  = 

This  split-operator  approach  is  sometimes  used  to  solve  stiff  systems  of  equations  (i.e. 
systems  with  widely  differing  time  scales)  such  as  those  that  result  when  hnite-rate 
chemistry  is  included  in  a  CFD  solution  (24,  25). 

An  approximate  numerical  solution  can  be  obtained  using  a  Strang  symmet¬ 
ric  splitting  approach  (26,  24,  23).  Given  difference  operators  AIat  and  £at  that 
are  second-order  accurate  for  Eqs.  2.22a  and  2.22b  respectively,  then  the  composed 
operator  Dat/2A^At'^^At/2  is  second-order  accurate  for  Eq.  2.15  (26).  Furthermore, 
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the  stability  characteristics  of  the  composed  operator  are  determined  by  the  stability 
characteristics  of  the  individual  operators  Af  At  and  £at  (26,  23).  By  choosing  ap¬ 
propriate  discretizations  for  Eqs.  2.22a  and  2.22b,  the  large- iV  time  step  restriction 
is  greatly  reduced  or  eliminated. 

To  take  advantage  of  the  efficiencies  of  the  multi-domain  approach,  Eq.  2.22a 
is  transformed  to  the  time  domain  in  a  manner  similar  to  Eq.  2.18.  Together  with 
Eq.  2.22b,  this  transformed  equation  becomes  the  split-domain  harmonic  balance 
form  of  Eq.  1.1. 


dj  dm  ^ 

dr  dx 

|  +  S(|)  =  0. 


(2.23a) 

(2.23b) 


The 
time  level 

step  1: 

step  2: 
step  3: 
step  4: 
step  5: 


steps  required  to  advance  the  solution  one  iteration  from  time  level  n  to 
n  +  1  are: 

Given  the  solution  vector  at  time  level  n,  advance  Eq.  2.23b  one-half 
pseudo-time  step  to  produce 

Calculate 

With  ^  ,  advance  Eq.  2.23a  one  full  pseudo-time  step  to  produce  ^ 
Calculate 

With  ,  advance  Eq.  2.23b  one-half  pseudo-time  step  to  produce 


The  split-domain  harmonic  balance  approach  has  several  advantages  over  pre¬ 
vious  harmonic  balance  implementations.  First  and  foremost,  it  greatly  reduces  the 
stability  restriction  for  large  N.  Secondly,  it  requires  only  one  DFT/IDFT  calcu¬ 
lation  per  point  per  iteration,  regardless  of  the  time  integration  scheme  employed. 
Another  advantage  lies  in  the  fact  that  the  system  of  equations  represented  by  the 
time-domain  PDE,  Eq.  2.23a,  is  uncoupled,  and  each  equation  has  the  same  form  as 
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the  original  conservation  equation.  Thus  a  solution  scheme  developed  for  the  orig¬ 
inal  equation  can  be  applied  in-turn,  without  modihcation,  to  each  equation  in  the 
harmonic  balance  system.  This  makes  the  split-domain  approach  easy  to  implement 
in  an  existing  code. 

The  split-domain  harmonic  balance  formulation  has  one  potential  disadvan¬ 
tage.  As  a  result  of  splitting  the  harmonic  balance  equation,  a  new  discretization 
error  is  introduced.  This  error  causes  the  steady-state  solution  to  be  dependent  on 
the  numerical  time-integration  step  size.  At  (25).  An  investigation  of  the  impact 
of  splitting  error  on  split-domain  harmonic  balance  solutions  was  conducted  (Ap¬ 
pendix  B)  and  it  was  found  that  splitting  error  is  not  a  signihcant  factor  for  the 
solver  implementations  studied  in  this  effort.  It  could  become  a  factor,  however,  for 
an  implementations  that  allow  very  large  numerical  time  steps. 

Like  the  time-domain  PDE,  the  system  of  equations  represented  by  the  fre¬ 
quency  domain  ODE,  Eq.  2.23b,  is  also  uncoupled,  with  each  equation  having  the 
form 


— i  fco;  Cfc  =  0,  0  <  k  <  N.  (2.24) 

dr 

Eq.  2.24  is  a  linear  ODE  with  an  easily  obtained  exact  solution.  Given  a  solution 
at  time  level  n  as  an  initial  condition,  the  exact  solution  to  Eq.  2.23b  at  pseudo-time 
Ar/2  is  given  by 


,  {)<k<N.  (2.25) 

Insight  into  the  split-domain  harmonic  balance  iterative  solution  process  is  gained 
by  recognizing  that  Eq.  2.25  represents  a  small  shift  of  the  solution  vector  elements 
in  physical  time.  The  time-shift  property  for  the  discrete  Fourier  transform  states 
that,  given  the  time  sampled  sequence  ^  with  length  2N  + 1,  and  its  discrete  Fourier 
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transform  then  (27) 

—  ^  _ i27Tna 

^k-a^Cne  2Ar+l, 

where  the  subscripts  k  and  n  range  over  the  elements  in  the  vector 
side  of  Eq.  2.25  is  similar  to  that  of  Eq.  2.26,  provided 

(2N  +  l)ci;Ar 
a  = - - - . 

Att 

Thus  the  combined  result  of  the  operations  represented  by  the  split-domain  solution 
steps  4,  5,  1,  and  2  (in  that  order)  approximates  a  small  physical-time  shift  of 
the  elements  in  the  vector  ^  along  the  continuous  periodic  function  given  by  the 
approximating  Fourier  series,  Eq.  2.2. 

With  this  understanding,  the  split-domain  harmonic  balance  solution  process 
can  be  interpreted  as  a  sequence  of  forward  integrations  of  Eq.  2.23a  in  pseudo-time 
(which,  because  Eq.  2.23a  has  the  same  form  as  Eq.  2.1,  is  just  a  scaled  integration 
in  physical  time),  surrounded  by  backward  physical-time  shifts  of  the  elements  in  the 
solution  vector  (Fig.  2.1).  The  shift-integrate-shift  sequence  is  continued  until  the 
end  values  are  the  same  as  the  starting  values,  at  which  point  the  overall  steady-state 
solution  is  reached. 

2.5  Adaptive  Split-Domain  Harmonic  Balance 

The  adaptive  split-domain  harmonic  balance  method  minimizes  the  computa¬ 
tional  cost  of  the  harmonic  balance  solution  by  automatically  tailoring  the  number 
of  Fourier  frequencies  included  in  the  solution  according  to  the  flow  at  a  given  point, 
on  a  point-by-point  (or  for  a  finite  volume  discretization,  cell-by-cell)  basis.  This  is 
accomplished  by  means  of  a  frequency  augmentation  approach.  With  this  approach, 
the  solution  is  begun  with  a  user-specihed  minimum  initial  number  of  frequencies. 
As  the  solution  develops,  frequencies  are  added  in  hxed  increments  to  individual  grid 
points  until  a  hnal  frequency  distribution  and  solution  are  obtained. 


(2.26) 

.  The  right  hand 

(2.27) 


2-12 


■  Space 

Figure  2.1  Time-sample  Shifting  that  Occurs  with  Each  Iteration  of  the  Split- 
domain  Solution  Process 

The  frequency  augmentation  approach  was  chosen  because  it  is  simple  and  com¬ 
putationally  efficient,  especially  for  problems  that  are  elliptic  in  time  (e.g.  subsonic 
flow  problems).  Computational  efficiency  stems  from  the  fact  that  the  low-frequency 
coefficients  included  in  the  early  solutions  tend  to  have  low-frequency  spatial  varia¬ 
tions  (Fig.  2.2),  and  thus  take  the  most  time  to  converge  to  steady  state.  By  begin¬ 
ning  with  a  small  number  of  Fourier  frequencies,  much  of  the  work  of  converging  the 
low-frequency  coefficients  is  accomplished  while  solving  a  reduced  problem.  Subse¬ 
quent  addition  of  higher  Fourier  frequencies  introduces  mostly  high-spatial- frequency 
errors  which  are  quickly  removed  from  the  solution.  In  this  respect,  frequency  aug¬ 
mentation  is  similar  to  a  full  multigrid  convergence  acceleration  method,  in  which 
a  solution  is  converged  on  a  series  of  computational  grids,  beginning  with  a  coarse 
grid  and  continuing  on  successively  hner  grids  until  a  hnal  solution  is  obtained. 

To  implement  a  frequency  augmentation  approach,  it  is  necessary  to  identify 
which  points  in  the  computational  grid  need  augmentation,  and  to  properly  handle 
the  frequency  transitions  that  occur  when  a  cell  and  its  neighbor  are  solved  with 
different  numbers  of  frequencies.  Determining  which  cells  need  additional  Fourier 
frequencies  requires  a  means  of  measuring  the  quality  of  the  solution,  and  some 
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Figure  2.2  Illustration  of  Typical  Spatial  Variation  of  a  Low-frequency  (top)  and 
High-frequency  (bottom)  Fourier  Coefficient 

criteria  against  which  that  quality  can  be  compared.  It  also  requires  a  strategy  that 
governs  how  often  the  criteria  are  applied.  Handling  frequency  transitions  requires 
the  development  of  a  robust  multi-frequency  approach  for  solving  the  time-domain 
equation,  Eq.  2.23a. 

2.5.1  Frequency  Augmentation  Criteria  and  Scheduling.  The  frequency 
augmentation  approach  requires  a  reliable  indicator  of  solution  fidelity  at  each  com¬ 
putational  cell.  Since  the  final  time-response  of  the  flow  is  not  known  a  priori^  the 
indicator  must  rely  only  on  the  current  (and  possibly  past)  state  of  the  solution. 
The  indicator  chosen  for  this  research  is  the  fraction  of  spectral  energy  contained  in 
the  highest  computed  Fourier  Frequency,  given  by 

En  =  Vf  (2.28) 

^n=0 

The  assumptions  underlying  the  choice  of  Etv  as  an  adaptation  metric  are  that  the 
majority  of  the  spectral  energy  is  contained  in  the  low-frequency  terms,  and  that  the 
fraction  of  energy  contained  in  the  highest  calculated  frequency  decreases  as  more 
terms  are  included  in  the  approximating  series.  Thus  the  quality  of  the  solution  can 
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be  gauged  by  the  amount  of  energy  in  the  highest  frequency  term.  Physically,  these 
assumptions  require  that,  aside  from  localized  discontinuities,  the  modeled  flow  held 
is  smoothly  varying  throughout  its  period.  In  addition,  if  discontinuities  do  appear, 
they  must  not  be  impulse  discontinuities,  i.e.,  the  discontinuity  should  be  a  step 
rather  than  a  spike.  These  physical  requirements  are  consistent  with  the  how  helds 
of  interest,  and  proved  valid  for  all  test  cases  examined. 

The  decision  to  augment  frequencies  at  a  point  is  made  by  comparing  to 
a  threshold  value,  Ethresh-  Because  Ej^  tends  to  mirror  the  spatially  oscillatory  na¬ 
ture  of  the  high-frequency  coefficients,  it  is  smoothed  with  an  unweighted  5-point 
spatial  average  prior  to  thresholding.  When  the  smoothed  Ejsi  exceeds  the  thresh¬ 
old,  additional  frequencies  are  incorporated  into  the  solution  at  that  point.  Fourier 
coefficients  for  the  new  frequencies  are  initialized  to  zero.  In  the  case  of  a  vector 
equation,  Ejsi  is  calculated  for  each  of  the  variables  expanded  in  a  series,  and  the 
solution  is  augmented  if  any  one  of  these  exceeds  Ethresh-  Selection  of  Ethresh  is  based 
on  experience  and  the  desired  solution  fidelity.  See  Sections  4.4  and  5.4  for  more  on 
the  impact  of  threshold  value  on  adapted  solutions. 

Threshold-based  augmentation  was  supplemented  by  two  forms  of  non-threshold 
based  augmentation.  The  hrst  of  these  was  fringe  augmentation.  The  purpose  of 
fringe  augmentation  was  to  increase  the  size  of  a  threshold-augmented  region.  In 
some  test  cases  with  very  rapid  transitions  between  smooth  and  discontinuous  flow, 
the  location  of  the  transitions  changed  as  frequencies  were  added  and  the  solution 
was  rehned.  In  those  situations,  it  was  sometimes  necessary  to  augment  a  small 
fringe  of  cells  adjacent  to  threshold-identihed  cells  in  order  to  allow  the  transition  to 
shift  in  the  direction  of  lower  frequency  content. 

The  second  non-threshold-based  augmentation  was  designed  to  enforce  a  min¬ 
imum  number  of  consecutive  cells  with  the  same  frequency  content.  In  some  cases, 
threshold-based  augmentation  can  result  in  small  (1-2  cell)  segments  of  the  compu¬ 
tational  domain  having  a  different  frequency  content  than  their  neighbors.  To  avoid 


2-15 


New  Frequency  Distribution. 


Figure  2.3  Effect  of  Fringe  Augmentation  and  Pixelation  on  a  Threshold-based 
Frequency  Distribution 

this  situation,  a  pixelation  process  was  applied  to  each  newly  augmented  frequency 
distribution.  In  this  process,  the  grid  was  divided  into  contiguous,  non-overlapping 
segments  with  a  uniform  size,  or  pixelation  width.  Within  each  segment,  the  fre¬ 
quency  content  was  set  to  the  maximum  found  within  that  segment.  This  guaranteed 
a  minimum  number  of  contiguous  cells  with  the  same  frequency  content,  while  also 
guaranteeing  frequency  content  greater  than  or  equal  to  that  required  to  meet  the 
augmentation  threshold.  The  cumulative  effect  of  both  fringe  augmentation  and 
pixelation  is  illustrated  in  Fig.  2.3. 

For  the  current  research,  both  the  augmentation  fringe  width  and  pixelation 
width  were  controlled  by  user  input  at  run  time.  Each  could  be  disabled  when  not 
needed  by  setting  the  fringe  width  to  zero,  and  the  pixelation  width  to  one. 

Once  a  cell  was  identified  for  augmentation,  its  frequency  content  was  increased 
by  a  predetermined  increment,  chosen  to  minimize  compute  time.  The  majority  of 
the  run  time  required  to  solve  the  split-domain  harmonic  balance  equations  is  com¬ 
posed  of  two  components-the  time  associated  with  solving  the  time-domain  equations 
(Eq.  2.23a),  and  the  time  spent  performing  the  necessary  Fourier  transforms.  Run 
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time  associated  with  the  time-domain  equations  increases  linearly  with  N.  There 
is  no  performance  advantage  favoring  any  particular  increment,  so  long  as  the  in¬ 
crement  is  relatively  small.  There  can,  however,  be  a  signihcant  difference  in  the 
run  time  of  the  Fourier  transforms  for  different  N,  as  discussed  in  Section  2.2.  Nu¬ 
merical  tests  were  conducted  to  determine  values  of  N  which  gave  the  best  overall 
performance.  The  results  of  these  tests  are  documented  in  Appendix  C. 

Like  the  augmentation  increment,  the  timing  of  frequency  adaptation  was 
based  on  run-time  performance  considerations.  Frequency  augmentation  is  most 
effective  when  it  is  based  on  solutions  that  are  representative  of  the  hnal  solution. 
If  adaptation  is  attempted  too  early  in  the  solution  process,  unneeded  frequencies 
could  be  added  based  on  transient  flow  structures  that  are  not  present  in  the  hnal 
solution.  On  the  other  hand,  if  the  solution  is  allowed  to  develop  too  long  before 
adaptation  is  attempted,  the  result  could  be  wasted  computational  effort.  The  goal 
of  adaptation  scheduling  is  to  identify  the  “right”  times  to  adapt  the  solution. 

A  dual-trigger  adaptation  scheduling  strategy  is  taken.  The  primary  trigger  is 
based  on  a  modihed  L2  norm  of  the  change  in  the  solution  during  one  iteration.  The 
L2  norm,  or  residual,  is  dehned  as 


R 


m 


(2.29) 


The  L2  norm  is  a  measure  of  the  remaining  error  in  the  solution.  In  Eq.  2.29,  ni 
is  the  number  of  cells  in  the  computational  grid,  A,^  is  the  change  in  the  solution 
vector  in  one  iteration,  and  the  overset  ~  indicates  complex  conjugation.  Adaptation 
is  triggered  when  log^Q(i?)  drops  by  a  user-specihed  amount,  indicating  that  a  desired 
level  of  solution  development  has  been  reached.  The  secondary  trigger  is  based  on  the 
number  of  iterations  completed.  The  iteration-based  scheduling  serves  as  a  backup 
to  residual-based  scheduling  in  the  rare  cases  where  solution  convergence  stalls  and 
the  specihed  residual  drop  is  not  achieved. 
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The  residual  drop  and  iteration  count  are  measured  relative  to  a  reference 
residual  and  iteration  number.  The  initial  reference  values  are  set  after  the  hrst 
iteration  of  the  solution.  The  reference  values  are  then  reset  whenever  the  solution 
is  adapted.  Adaptation  trigger  values  for  the  first  and  subsequent  adaptations  are 
different.  The  initial  trigger  values  are  set  to  allow  time  for  the  solution  to  develop 
from  a  poor  initial  guess.  Subsequent  trigger  values  are  set  to  allow  errors  introduced 
by  frequency  augmentation  to  be  removed  from  the  solution,  and  to  further  rehne 
the  solution  a  small  amount.  The  initial  trigger  values  are  typically  much  larger  than 
the  subsequent  trigger  values. 

The  adaptive  frequency  augmentation  algorithm  is  summarized  in  Fig.  2.4. 

2.5.2  Treatment  of  Frequeney  Transitions.  Discretization  of  the  spatial 
derivative  in  the  time-domain  portion  of  the  split-domain  harmonic  balance  equa¬ 
tions  (Eq.  2.23a)  requires  the  addition  and/or  subtraction  of  the  solution  vector  in 
a  cell,  (or  the  corresponding  fluxes,  with  the  other  solution  vectors  in 

its  discretization  stencil.  This  presents  a  problem  when  those  cells  have  different 
maximum  frequencies.  Not  only  do  the  solution  vectors  have  different  numbers  of 
elements,  but  those  elements  correspond  to  the  state  of  the  flow  at  different  points 
in  time;  they  have  different  sample  rates.  To  solve  this  problem,  the  solution  vectors 
must  be  resampled  so  that  the  sample  rate  is  consistent  across  the  discretization 
stencil  (Fig.  2.5). 

Two  different  resampling  methods  are  employed  in  the  present  implementa¬ 
tion.  These  methods  include  truncation/zero-padding  in  the  frequency  domain, 
and  linear  interpolation  in  the  time  domain.  The  primary  means  of  resampling 
is  truncation/zero-padding  in  the  frequency  domain.  With  this  method,  a  vector  is 
down-sampled  by  taking  its  Fourier  transform,  truncating  the  results,  and  transform¬ 
ing  the  truncated  Fourier  coefficient  vector  back  to  the  time  domain.  Upsampling  is 
achieved  in  a  similar  manner,  except  that  instead  of  truncating  the  Fourier  coefficient 
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Figure  2.4 


Flowchart  for  the  Adaptive  Harmonic  Balance  Algorithm 
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Figure  2.5  Sample-rate  Interpolation  Required  at  a  Frequency  Transition  from 
TV  =  2  to  iV  =  3 


vector,  it  is  padded  with  new  zero-valued  high-frequency  coefficients.  This  approach 
results  in  very  smooth  interpolation.  When  the  results  are  used  in  a  linear  operation 
such  as  addition  or  subtraction,  the  approach  is  equivalent  to  performing  the  same 
linear  operation  directly  in  the  frequency  domain. 

One  drawback  with  frequency-domain  truncation/zero  padding  is  that  the  in¬ 
terpolated  values  are  not  bounded  by  the  original  data.  When  sample  rates  are  very 
small,  interpolation  can  result  in  non-physical  values,  as  demonstrated  in  Fig.  2.6. 
When  this  occurs,  linear  interpolation  in  the  time-domain  is  applied. 

With  the  linear  interpolation  method,  conservative  variables  are  obtained  at 
the  required  sample  times  by  linearly  interpolating  the  calculated  values.  Besides 
guaranteeing  that  the  interpolated  values  are  bounded  by  the  computed  values,  this 
approach  has  the  advantage  that  it  can  be  applied  entirely  in  the  time  domain,  and 
does  not  require  any  additional  Fourier  transforms. 

Despite  the  higher  computational  cost  and  potential  low-sample-rate  problems, 
frequency-domain  truncation/zero-padding  was  selected  as  the  primary  resampling 


2-20 


0  0.2  0.4  0.6  0.8  1 

Normalized  Time 


Figure  2.6  Non-physical  Interpolated  Values  Produced  by  Frequency-domain  Zero¬ 
padding 

method.  This  choice  was  based  on  the  quality  of  the  solution  computed  at  a  sample 
rate  transition.  The  superiority  of  the  method  over  linear  interpolation  is  illustrated 
in  Fig.  2.7.  This  hgure  contains  plots  of  contours  of  constant  Fourier  coefficient 
magnitude  for  both  types  of  interpolation  at  a  transition  from  7  frequencies  (15 
samples/period)  to  16  frequencies  (33  samples/period).  Also  included  in  the  plots 
(dashed)  are  coefficient  magnitude  contours  for  a  solution  with  a  constant  16  frequen¬ 
cies.  The  frequency-domain  truncation/zero  padded  solution  is  much  smoother  and 
more  closely  matches  the  constant-sample-rate  results  than  the  linearly-interpolated 
solution. 

In  the  adaptive  split-domain  harmonic  balance  solver  written  for  this  research, 
a  3-stage  Runge-Kutta  time-integration  scheme  was  used  to  advance  the  time-domain 
PDF,  Eq.  2.23a  (See  Section  4.2).  In  order  to  achieve  a  smooth  solution  at  frequency 
transitions,  it  was  necessary  to  resample  the  transition  boundaries  at  every  integra¬ 
tion  stage.  Attempts  to  time-lag  the  transitions  by  freezing  the  resampled  values 
at  the  hrst  integration  stage  resulted  in  discontinuities  at  the  transition  point,  even 
when  both  sides  of  the  transition  had  the  same  number  of  frequencies  and  no  inter¬ 
polation  was  required. 
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Linear  Interploation  Frequency  T  runcation/Zero-Padding 


Figure  2.7  Comparison  of  Fourier  Coefficient  Magnitudes  at  a  Transition 
from  7  Frequencies  to  16  Frequencies,  for  Solutions  Computed 
with  Time-domain  Linear  Interpolation  (left)  and  Frequency- domain 
Truncation/Zero-padding  (right)  Resampling  Methods 


As  a  consequence  of  the  need  to  resample  the  transition  boundaries  at  each  inte¬ 
gration  stage,  the  previously  uncoupled  system  of  equations  represented  by  Eq.  2.23a 
becomes  effectively  coupled  through  the  resampling  operation.  This  signihcantly  in¬ 
creases  the  storage  requirements  of  the  solver.  When  the  time-domain  system  of 
equations  was  uncoupled,  the  numerical  integration  scheme  could  be  applied  sepa¬ 
rately  to  each  equation  in  the  system,  reducing  the  storage  requirement  for  inter¬ 
mediate  solutions  to  a  single  time  sample.  But  when  the  equations  are  coupled, 
intermediate  solutions  for  all  samples  must  be  stored.  For  this  reason,  it  is  impor¬ 
tant  that  the  adaptive  split-domain  harmonic  balance  method  be  implemented  with 
a  numerical  integration  scheme  with  low  storage  requirements  if  that  scheme  is  a 
mulit-step  scheme. 


2.6  Multigrid  and  Adaptive  Split-Domain  Harmonic  Balance 

The  Full  Approximation  Storage  (FAS)  multigrid  method  is  a  commonly  used 
technique  for  accelerating  the  calculation  of  steady-state  solutions  to  nonlinear  con¬ 
servation  equations.  The  method  works  by  transferring  a  partially  converged  solution 
to  a  coarser  computational  grid,  where  a  coarse  grid  correction  is  calculated.  This 
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correction  is  then  applied  to  the  solution  on  the  fine  grid.  Because  the  coarse  grid 
contains  fewer  grid  points,  calculation  of  the  coarse  grid  correction  is  faster  than 
calculation  of  the  solution  on  the  fine  grid. 

Implementation  of  FAS  multigrid  in  an  adaptive  split-domain  harmonic  balance 
solver  requires  some  special  considerations.  In  the  following  discussion,  the  basic 
FAS  multigrid  scheme  as  found  in  (28)  is  presented,  followed  by  the  specifics  of  the 
adaptive  split-domain  implementation. 

Theory  of  FAS  Multigrid.  Consider  an  equation  of  the  form 

j^hjjh  ^  ^2.30) 

defined  on  a  grid  with  spacing  h,  where  is  a  nonlinear  operator,  is  the  exact 
solution,  and  is  a  forcing  function.  Let  be  an  approximate  solution  to  Eq.  2.30, 
and  be  the  error  in  the  approximate  solution.  The  goal  of  the  FAS 

multigrid  scheme  is  to  quickly  calculate  an  estimate  of  the  solution  error, 

Substituting  +  u^  into  Eq.  2.30  and  subtracting  from  both  sides 

gives 

L^{V^ +  u^)- -  L^u^.  (2.31) 

If  the  terms  in  Eq.  2.31  are  sufficiently  smooth,  the  equation  can  be  transferred  or 
“restricted”  onto  a  coarse  grid  with  grid  spacing  2h  without  much  loss  of  accuracy. 
Given  a  coarse  grid  operator,  and  a  restriction  operator,  the  restricted 
equation  is  given  by 

l2h^l2h^h  ^  y2h^  _  l}h^j2h^h^  ^  J2h^j^h  _  j^h^hy  ^2.32) 

The  coarse  grid  operator  may  or  may  not  be  the  same  as  the  fine  grid  operator, 
L^.  is  the  approximate  solution  error  on  the  coarse  grid.  Eq.  2.32  can  be 
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rewritten  to  have  the  same  form  as  Eq.  2.30,  i.e., 


L 


2h^2h 


f 


2h 


(2.33) 


where 


(2.34) 

j2h  ^  j2h^fh  _  ^  L^^(ll^u^).  (2.35) 

The  FAS  multigrid  scheme  solves  for  on  the  coarse  grid,  which,  because  it  has 
half  as  many  points  as  the  hne  grid,  requires  half  as  much  work.  In  addition,  since 
the  coarse  grid  spacing  is  larger,  in  many  cases  a  larger  integration  time  step  may  be 
taken.  Once  a  solution  for  is  obtained,  the  error  in  the  solution,  is  recovered 
and  interpolated  to  the  hne  grid  to  form  a  coarse  grid  correction.  The  hne  grid 
solution  is  updated  by  adding  the  coarse  grid  correction  to  the  original  approximate 
solution,  as  illustrated  in  Eq.  2.36. 

(«")„ew  =  («")old  +  <d)  (2.36) 

Because  Eq.  2.33  has  the  same  form  as  Eq.  2.30,  the  FAS  multigrid  algorithm 
is  easily  extended  to  more  than  two  levels.  Each  coarse  grid  becomes  the  hne  grid 
for  a  still  coarser  grid.  In  theory,  this  can  continue  until  there  is  only  one  point  in 
the  interior  of  the  coarsest  grid. 

A  typical  multigrid  implementation  begins  by  performing  a  small  number  of 
iterations  on  the  hnest  grid  to  remove  high  spatial-frequency  errors.  This  smoothed 
solution  is  restricted  to  the  next  coarsest  grid,  where  the  solution  is  again  smoothed 
for  a  small  number  of  iterations.  The  process  of  smooth  and  restrict  is  repeated 
until  the  coarsest  grid  is  reached,  at  which  point  the  solution  is  converged  to  steady 
state.  Once  steady  state  is  reached,  a  coarse  grid  correction  is  made  to  the  next  hner 
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grid,  and  a  small  number  of  iterations  are  performed  to  remove  errors  introduced  by 
prolongation.  The  process  is  repeated  until  the  hnest  grid  is  reached,  completing  one 
multigrid  cycle.  This  down-and-up  pattern  is  referred  to  as  a  “V”  cycle  (Fig.  2.8). 
Other  cycles  with  different  restriction/prolongation  patterns  (e.g.,  the  “W”  cycle) 
are  also  sometimes  used. 


Figure  2.8  Illustration  of  a  4-Level  Multigrid  V  Cycle 


FAS  Multigrid  and  Split-domain  Harmonic  Balance.  In  order  to  apply  FAS 
Multigrid  to  the  split-domain  harmonic  balance  equations,  they  must  be  put  in  the 
form  of  Eq.  2.30.  Let  Lj  be  dehned  as  the  time-integration  operator  that  advances  the 
solution  to  the  frequency- domain  ODE  (Eq.  2.23b)  by  one  half  time  step.  Similarly, 
let  Lt  be  defined  as  the  time-integration  operator  that  advances  the  time-domain 
PDE  by  one  full  time  step.  The  process  of  integrating  the  solution  vector  ^  from 
time  level  n  to  time  level  n  -|-  1  with  the  split-domain  scheme  can  then  be  written  as 

^^^+1  =  (2.37) 

Subtracting  from  both  sides  and  dividing  by  the  integration  time  step.  At,  gives 

Pn+l  _  Cn  1 

=  -^ILfTLT-'Lf  -  /|r,  (2.38) 
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where  /  is  the  identity  operator.  The  left-hand  side  of  Eq.  2.38  is  a  hrst-order 
forward-difference  approximation  of  the  pseudo-time  derivative  ^  added  to  the 
steady-state  harmonic  balance  equation  to  allow  the  use  of  a  time-marching  solution 
method.  Removing  the  derivative  results  in  an  operator-notation  expression  for  the 
steady-state  split-domain  harmonic  balance  equations  that  is  in  the  desired  form, 
i.e.  =  f,  where  /  =  0,  and 

(2.39) 

To  calculate  a  coarse  grid  correction  for  the  split-domain  harmonic  balance 
equations  at  time  level  n,  the  updated  solution  at  time  level  n  -|-  1,  must  first 
be  calculated.  The  change  in  the  solution  vector,  is  computed  and  divided 

by  the  pseudo-time  step.  At.  The  result,  along  with  the  solution  at  time  level  n,  is 
restricted  to  the  next  coarsest  grid  and  the  coarse  grid  correction  is  calculated.  The 
correction  is  added  to  The  previously  calculated  is  discarded. 

FAS  Multigrid  and  Frequency  Augmentation.  To  implement  adaptive  frequency 
augmentation  and  FAS  multigrid  together,  no  major  changes  are  required  in  either 
method.  A  strategy  is  required  to  govern  what  grid  levels  augmentation  will  occur 
on,  and  how  frequency  maps  will  be  transfered  from  one  grid  level  to  another.  The 
only  additional  requirement  is  that  frequency  transitions  must  be  properly  handled 
during  restriction  and  prolongation. 

In  the  current  research,  frequency  augmentation  was  applied  only  on  the  hnest 
grid.  Frequency  map  consistency  from  the  hnest  grid  to  the  coarsest  grid  was  main¬ 
tained  by  setting  the  frequency  map  pixelation  width  to  2^“^,  where  k  is  the  depth  of 
the  multigrid  cycle.  This  width  guarantees  that  the  two  hne  grid  cells  contributing 
to  a  coarse  grid  cell  both  have  the  same  sample  rate  for  all  grid  levels  (Fig.  2.9). 
Thus  every  region  of  the  grid  is  represented  by  the  same  frequency  content  on  all 
grid  levels.  This  approach  eliminates  frequency  transitions  during  restriction,  be- 
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Figure  2.9  Illustration  of  How  a  Pixelation  Width  of  2^“^  Preserves  Spatial  Fre¬ 
quency  Content  and  Eliminates  Frequency  Transitions  During  Restric¬ 
tion 

cause  all  frequency  transitions  occur  between  cells  that  contribute  to  different  coarse 
grid  cells.  Frequency  transitions  are  still  encountered  during  prolongation,  however. 
These  transitions  are  handled  by  zero-padding  the  high-frequency  Fourier  coefficients 
of  the  cell  with  the  smaller  sample  rate. 


2-27 


III.  Application  of  Split-Domain  Harmonic  Balance  to  Burgers  ’ 

Equation 

3. 1  Introduction 

The  purpose  of  the  research  documented  in  this  chapter  was  to  examine  the 
behavior  of  the  split-domain  harmonic  balance  method  when  applied  to  problems 
with  strong  discontinuities  moving  over  large  regions  of  the  grid.  A  useful  model 
for  this  research  was  the  1-D  inviscid  Burgers’  equation,  a  simplihed  form  of  Euler’s 
equation  that  yields  traveling  discontinuities  in  the  flow  held  for  large  amplitude 
periodic  disturbance  boundary  conditions.  The  study  included  the  application  of 
the  harmonic  balance  method  to  two  families  of  unsteady  boundary  conditions  -  one 
based  on  a  single-frequency  sine  wave  and  the  second  based  on  a  simulated  wake 
function.  The  amplitude  and  frequency  of  the  boundary  conditions  were  varied  to 
generate  test  cases  with  a  wide  range  of  how  properties,  from  smooth  and  continuous 
to  strongly  discontinuous.  The  ehect  of  number  of  Fourier  frequencies  included  in 
the  harmonic  balance  calculation  of  the  solution  to  these  diherent  test  cases  was 
of  particular  interest.  To  provide  a  basis  for  comparison  of  accuracy,  stability,  and 
performance,  harmonic  balance  solvers  based  on  prior  harmonic  balance  approaches 
were  also  implemented  and  tested. 


3.2  Solver  Implementation 

Four  harmonic  balance  CFD  solvers  were  written  to  solve  the  one-dimensional 
inviscid  Burgers  equation,  given  by 


du  1  dE 
dt  2  dx 


(3.1) 


where  u{x,t)  is  the  dependent  variable,  t  is  the  temporal  variable,  and  x  is  the  spa¬ 
tial  variable.  The  four  solvers  included  a  direct  harmonic  balance  implementation 


3-1 


(Eq.  2.15),  a  time-domain  implementation  (Eq.  2.18),  a  frequency-domain  imple¬ 
mentation  (Eq.  2.21),  and  a  split-domain  implementation  (Eq.  2.23). 

The  direct,  time-domain,  and  frequency-domain  solvers  were  implemented  with 
an  explicit  MacCormack  discretization  scheme,  modihed  to  incorporate  source  terms 
(24).  MacCormack’s  scheme  is  a  two-step  solution  scheme  that  is  second  order  in 
both  time  and  space.  For  a  scalar  conservation  equation  with  source  term,  the 
modihed  MacCormack  scheme  is  given  by 


cn+l 

>2 


C  -  Ai 


4(S‘+i)  -  «(«.“) 
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+  s(a 


|>(C)  -  <i>(g-i) 
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+  s(e; 


(3.2a) 


(3.2b) 


Discretizations  of  the  direct,  time-domain,  and  frequency- domain  harmonic  balance 
equations  were  obtained  by  substitution  of  the  appropriate  solution  vectors,  hux 
functions,  and  source  terms  into  Eqs.  3.2a  and  3.2b.  (See  Appendix  A  for  a  derivation 
of  the  harmonic  balance  Burgers’  hux  function  and  source  term  required  for  the  direct 
solver.)  For  the  frequency-domain  formulation,  substitution  resulted  in  Eqs.  3.3a  and 
3.3b,  and  for  the  time-domain  formulation,  Eqs.  3.4a  and  3.4b. 
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In  these  expressions,  the  squaring  of  a  vector  quantity  denotes  element-wise  squaring. 
The  quantities  and  S  are  as  dehned  in  Section  2.2.  The  residual  function  in  the 
time-domain  harmonic  balance  equation  (Eq.  2.21)  has  a  different  dehnition  for  each 
of  the  solution  stages,  and  is  dehned  as  the  expressions  in  square  brackets  in  Eqs.  3.4a 
and  3.4b.  Note  that  because  these  expressions  contain  only  linear  operations  (vector 
subtraction  and  scalar  division),  the  frequency  domain  solver  and  time-domain  solver 
should  be  equivalent. 

The  split-domain  solver  was  implemented  with  different  discretization  schemes 
for  the  time-domain  PDE,  Eq.  2.23a,  and  the  frequency-domain  ODE,  Eq.  2.23b. 
The  time-domain  PDE  was  solved  with  an  unmodihed  MacCormack  discretization, 
given  by 
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The  frequency-domain  ODE  was  solved  with  a  three-stage  Runge-Kutta  numerical 
integrator 
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(3.6) 


Note  that  Eq.  3.6  advances  the  ODE  by  one-half  pseudo-time  step,  as  required  by 
the  symmetric  Strang  splitting  approach  (see  Section  2.4). 

All  solvers  required  artihcial  dissipation  to  prevent  oscillations  near  discontinu¬ 
ities  in  the  solution.  Dissipation  was  incorporated  in  the  direct,  frequency-domain, 
and  split-domain  formulations  by  adding  second  derivative  smoothing  (Eq.  3.7)  of 
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the  Fourier  coefficients  at  the  end  of  each  iteration,  according  to 


Ci,  new 


—  +  Ci+l 

A.x^ 


(3.7) 


Dissipation  was  added  to  the  time-domain  solver  by  applying  the  same  smoothing 
operator  to  the  time-domain  solution  vector,  Due  to  the  linearity  of  the  Fourier 
transform  and  artihcial  dissipation  operators,  this  added  dissipation  was  mathemat¬ 
ically  equivalent  to  that  of  the  other  solvers.  The  dissipation  parameter  a  was  a 
small  constant,  of  order  l.Oe-6,  that  controlled  the  amount  of  applied  damping. 

Burgers’  equation  is  hyperbolic  in  time.  For  positive  values  of  the  dependent 
variable,  u,  all  flow  information  travels  in  the  direction  of  increasing  spatial  coordi¬ 
nate.  Boundary  conditions  are  therefore  implemented  by  hxing  the  inflow  boundary 
values,  and  extrapolating  the  outflow  boundary  values  from  upstream.  A  0**^  order 
extrapolation  was  implemented.  For  the  direct,  frequency-domain,  and  split-domain 
solvers,  boundary  conditions  were  applied  in  the  frequency  domain.  Boundary  con¬ 
ditions  were  applied  in  the  time  domain  for  the  time-domain  solver. 

To  accelerate  convergence  to  steady  state,  local  time  stepping  was  employed. 
Using  the  dehnition  of  the  Courant-Friedricks-Lewy  (CFL)  stability  limit  for  the 
MacCormack  discretization  (29),  the  maximum  time  step  at  each  point  was  deter¬ 
mined  according  to 

At* 

At,  =  CFL - ^  (3.8) 

max(^J 

where  CFL  was  a  user-specihed  value  less  than  or  equal  to  1.0,  and  max(,^J  was  the 
maximum  element  in  the  vector  of  time-sampled  dependent  variables  at  the  grid 
point. 


3.3  Test  Configuration 

Solutions  were  obtained  for  two  families  of  periodic  inflow  boundary  condi¬ 
tions.  The  hrst  family  was  constructed  by  varying  the  amplitude  and  frequency  of  a 


3-4 


sinusoidal  oscillation  about  a  mean  value  of  1.0,  i.e., 


m(0,  t)  =  1.0  +  asin(ct;t) 


(3.9) 


where  a  is  the  amplitude  of  the  oscillation,  and  oo  is  the  frequency  of  the  oscillation 
in  radians  per  second.  The  second  family  consisted  of  amplitude  and  frequency 
variations  of  a  simulated  wake  function,  dehned  by 


m(0, t) 


1.0  0<t<|^ 

—  2uj 

1.0  +  asin2(2o;t)  ^<t<^ 


The  basic  forms  of  the  input  families  are  shown  in  Fig.  3.1. 
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Figure  3.1  Sinusoidal  (left)  and  Wake  Function  (right)  Boundary  Condition  Wave¬ 
forms 

For  each  of  these  families,  combinations  of  three  amplitudes  and  three  frequen¬ 
cies  were  applied,  for  a  total  of  nine  inflow  conditions.  The  three  amplitudes  were 
a  =  0.1,  0.3,  and  0.5,  while  the  frequencies,  defined  in  terms  of  /  =  uj/2n,  were 
/  =  0.75,  1.5,  and  3.0.  Over  a  2-unit  grid,  these  amplitudes  produce  a  wide  range 
of  flow  behavior,  from  smooth  and  continuous  to  strongly  shocked.  The  frequencies 
chosen  are  typical  of  the  reduced  frequencies  found  in  axial  compressor  simulations 
(5,  6,  9). 
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Inflow  boundary  conditions  for  the  time-domain  solver  were  obtained  by  eval¬ 


uating  Eqs.  3.9  and  3.10  at  2N  +  1  points  in  time,  given  by 


27171 

uj{2N  +  l) 


0<n<2N.  (3.11) 


Inflow  boundary  conditions  for  the  direct,  frequency- domain,  and  split-domain  solvers 
were  obtained  by  computing  Fourier  series  coefficients  for  N  frequencies.  For  Fq.  3.9, 
the  real  Fourier  series  coefficients  are  given  by 


1.0 

(3.12a) 

0.0 

1  <  n  <  iV 

(3.12b) 

L 

n  =  1 

(3.12c) 

[o.o 

1  <  n  <  iV. 

The  real  Fourier  series  coefficients  for  Fq.  3.10  are  given  by 


Oq  — 


Ctn, 


bn  = 
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(3.13b) 

(3.13c) 


Complex  Fourier  series  coefficients  for  the  positive  frequencies  were  obtained 
from  the  real  Fourier  series  coefficient  by  means  of  the  following  relations 

Co  =  ^ao  (3.14a) 

Cn  =  ^{ttn-ibn)-  (3.14b) 


3-6 


3-4  Results 

3.4-1  Accuracy.  From  the  standpoint  of  accuracy,  the  four  harmonic  bal¬ 
ance  implementations  were  equivalent,  producing  solutions  that  were  effectively  iden¬ 
tical.  Representative  solutions  are  shown  for  each  input  variation  in  Figs.  3.2  and 
3.3.  These  hgures  compare  the  solution  at  t  =  0  relative  to  the  input  period  with  an 
equivalent  fully  developed  time-accurate  calculation.  The  harmonic  balance  solutions 
were  generated  with  48  Fourier  frequencies  on  a  501  point  grid.  The  time-accurate 
solutions  were  obtained  on  the  same  grid  using  a  validated  MacCormack  scheme  with 
the  same  artihcial  dissipation  used  for  the  harmonic  balance  equation,  Eq.  3.7. 

As  can  be  seen  in  Figs.  3.2  and  3.3,  the  input  boundary  conditions  resulted 
in  solutions  ranging  from  smooth  traveling  waves  to  strong  moving  discontinuities. 
In  all  but  two  cases,  the  48  frequency  harmonic  balance  solutions  were  comparable 
to  the  time-accurate  solutions.  The  two  exceptions  were  the  sine  input  cases  with 
amplitudes  a  =  0.5  and  0.3  at  the  lowest  disturbance  frequency.  These  cases  contain 
signihcant  high-frequency  oscillations  in  their  solutions.  It  will  be  shown  that  these 
cases  require  additional  terms  in  the  approximating  series. 

Effect  of  Series  Length.  To  determine  the  effect  of  series  length  on  the  accuracy 
of  the  harmonic  balance  method,  solutions  were  generated  for  each  input  condition 
with  series  lengths  ranging  from  2  to  48  frequencies.  A  quantitative  measure  of  the 
difference  between  each  of  these  solutions  and  an  equivalent  time-accurate  solution 
was  obtained  by  calculating  the  difference  root  mean  square  (RMS)  for  10  equally 
spaced  temporal  samples  spanning  one  period  of  the  disturbance.  These  10  RMS 
differences  were  averaged  to  obtain  a  solution  difference.  The  results  are  plotted  in 
Figs.  3.4  and  3.5. 

The  results  show  that  the  differences  between  the  harmonic  balance  and  time- 
accurate  solutions  did  not  go  to  zero,  but  were  asymptotic  with  respect  to  approx¬ 
imating  series  length.  In  each  case  where  a  good  solution  was  obtained,  there  was 
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Disturbance  Frequency 


Figure  3.2 


Comparison  of  Time-accurate  (Solid  Line)  and  48-Frequency  Harmonic 
Balance  (Symbols)  Solutions  for  the  Sine  Input  for  t  =  0.  Inset  plots 
Show  Dependent  Variable  Magnitude  vs.  Nondimensional  Distance 
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Disturbance  Frequency 


Figure  3.3  Comparison  of  Time-accurate  (Solid  line)  and  48-Frequency  Harmonic 
Balance(Symbols)  Solutions  for  the  Wake  Function  Input  for  t  =  0. 
Inset  Plots  Show  Dependent  Variable  Magnitude  vs.  Nondimensional 
Distance 
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10° 


Figure  3.4  RMS  Error  Between  Time-accurate  and  Harmonic  Balance  Solutions — 
Sine  Input 


Figure  3.5  RMS  Error  Between  Time-accurate  and  Harmonic  Balance  Solutions — 
Wake  Function  Input 
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Figure  3.6  Comparison  of  48-Frequency  and  97-Frequency  Solutions  for  a  =  0.5, 
/  =  0.75,  on  a  501  Point  Grid 

a  series  length  corresponding  to  some  Fourier  frequency  beyond  which  no  improve¬ 
ment  occurred.  This  frequency  is  hereafter  referred  to  as  the  asymptotic  frequency, 
and  the  associated  solution  is  called  the  asymptotic  solution.  Solutions  based  on 
fewer  frequencies  are  called  sub-asymptotic,  while  those  with  more  frequencies  are 
called  super-asymptotic.  It  is  clear  from  Fig.  3.4  that  the  two  cases  with  oscillatory 
solutions  are  sub-asymptotic,  and  thus  required  additional  frequencies  to  minimize 
error.  This  was  conhrmed  by  generating  a  97  frequency  solution  to  the  a  =  0.5, 
/  =  0.75  case,  which  is  compared  with  the  48  frequency  solution  in  Fig.  3.6.  In  the 
higher-frequency  solution,  most  of  the  high  frequency  oscillations  have  been  removed. 

The  fact  that  the  differences  between  the  harmonic  balance  and  time-accurate 
solutions  were  asymptotic  with  respect  to  series  length  does  not  mean  that  the 
harmonic  balance  solutions  did  not  continue  to  converge.  The  truncation  error  in 
the  harmonic  balance  solution  simply  became  insignihcant  compared  to  differences 
caused  by  other  factors  such  as  a  slight  difference  in  shock  location. 
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Both  disturbance  amplitude  and  disturbance  frequency  inflnenced  the  asymp¬ 
totic  freqnency.  As  amplitnde  increased,  the  asymptotic  freqnency  also  increased. 
This  was  attribnted  to  the  presence  of  stronger  discontinnities  in  the  larger  ampli- 
tnde  solntions.  In  contrast,  as  the  distnrbance  frequency  increased,  the  asymptotic 
freqnency  usually  decreased.  This  behavior  is  explained  qnalitatively. 

If  the  shape  and  period  of  the  time  response  of  Bnrgers’  eqnation  at  every 
point  in  the  grid  is  hxed,  the  response  then  has  a  hxed  frequency  content  with 
signihcant  energy  np  to  some  freqnency  oJmax-  In  the  harmonic  balance  solntion, 
^max  is  expressed  as  some  mnltiple  of  the  fnndamental  freqnency,  say  rimax^-  Then 
nmax  =  ^max/^  IS  inversely  proportional  to  u,  increasing  as  u  decreases. 

True  inverse  proportionality  was  observed  in  some  test  cases  (e.g.,  the  sine 
inpnt  with  freqnency  change  from  /  =  1.5  to  /  =  3.0,  all  amplitndes),  bnt  for 
the  majority  of  the  cases  the  asymptotic  frequency  was  less  than  that  predicted  by 
the  simple  model.  This  is  consistent  with  the  fact  that  a  discontinuity  in  the  time 
response  is  generally  sharper  for  higher  disturbance  freqnencies,  and  thus  Umax  is  not 
constant. 


Effect  of  Grid  Density.  All  of  the  harmonic  balance  solntions  became 
dissipative  to  some  degree  as  the  compntational  grid  was  coarsened.  The  impact  on 
solntion  qnality  depended  on  the  natnre  of  the  flow  held.  For  solutions  with  strong 
discontinnities  the  effect  was  relatively  minor,  and  sometimes  benehcial,  while  for 
smooth,  small-amplitude  solutions  the  effect  introdnced  severe  damping. 

One  example  of  a  benehcial  dissipative  ehect  is  shown  in  Fig.  3.7.  This  hgure 
shows  that  while  the  coarse  grid  solntion  experiences  some  smearing  of  the  shock, 
there  is  almost  complete  elimination  of  the  non-physical  oscillations  present  in  the 
49-freqnency  hne  grid  solution.  Eliminating  these  oscillations  on  a  hne  grid  wonld 
reqnire  the  use  of  a  much  longer  approximating  Fonrier  series. 
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Figure  3.7  Smoothing  Effect  of  Coarse  Grid  on  Harmonic  Balance  Solution 

An  example  of  unfavorable  effect  of  a  coarse  grid  is  illustrated  in  Fig.  3.8.  This 
hgure  compares  the  501-point  time-accurate  solution,  the  101-point  time-accurate  so¬ 
lution,  and  the  101  point  harmonic  balance  solution.  The  coarse-grid,  time-accurate 
solution  shows  some  degradation,  primarily  in  the  form  of  a  phase-lag  in  the  peaks  of 
the  solution.  In  contrast,  the  coarse  grid  damping  effect  caused  considerable  degra¬ 
dation  in  the  harmonic  balance  solution.  The  harmonic  balance  method  was  more 
sensitive  to  grid  density  than  the  time-accurate  method. 

The  effect  of  grid  density  in  both  cases  was  traced  to  the  Fourier  coefficients 
corresponding  to  the  higher  computed  frequencies.  Figure  3.9(a)  shows  the  variation 
in  magnitude  of  one  high  frequency  (n  =  47)  coefficient  for  the  a  =  0.5,  /  =  0.75  sine 
input  on  the  501  point  grid.  The  computed  coefficient  shows  rapid  oscillation  in  the 
spatial  dimension.  Figure  3.9(b)  shows  the  same  coefficient  calculated  on  the  101 
point  grid.  In  this  case,  the  coarse  grid  did  not  contain  sufficient  spatial  resolution 
to  capture  the  oscillations  in  magnitude,  and  the  magnitude  of  the  coefficient  was 
under-predicted.  The  impact  of  poorly  resolved  high-frequency  coefficients  depended 
on  the  relative  importance  of  those  frequencies  in  the  harmonic  balance  solution. 
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Figure  3.8  Severe  Damping  Effect  of  Coarse  Grid  on  Harmonic  Balance  Solution 

In  the  case  shown  in  Fig.  3.7,  most  of  the  energy  in  the  solution  was  contained 
in  relatively  low  frequencies,  and  coarse-grid  damping  of  the  highest  frequencies 
resulted  in  benehcial  smoothing.  For  the  case  shown  in  Fig.  3.8,  however,  the  damped 
frequencies  comprised  a  signihcant  part  of  the  solution,  and  the  overall  accuracy  was 
degraded.  The  generalized  loss  of  higher-frequency  information  caused  the  harmonic 
balance  method  to  be  more  sensitive  to  grid  density  than  the  time-accurate  method. 

The  results  suggest  that  while  grid  density  is  important  for  harmonic  balance 
solutions,  grid  density  is  even  more  important  for  smooth  solutions  that  require  fewer 
terms  in  the  approximating  series.  Problems  that  require  many  terms  may  require 
less  grid  resolution,  partially  offsetting  the  cost  of  the  additional  terms. 


3.4-2  Stability.  One  of  the  primary  motivations  for  the  split-domain  ap¬ 
proach  was  improved  high-frequency  stability.  The  superiority  of  the  split-domain 
approach  was  clearly  evident  in  the  test  results.  For  low-frequency  {N  <  10)  so¬ 
lutions,  all  of  the  harmonic  balance  solvers  exhibited  good  stability,  successfully 
computing  solutions  with  CFL  >  0.95.  As  the  number  of  frequencies  increased,  the 


3-14 


Figure  3.9  Comparison  of  Spatial  Variation  of  High-frequency  Coefficient  on  (a) 
501  Point  Grid  and  (b)  101  Point  Grid 

maximum  stable  CFL  of  the  direct,  time-domain,  and  frequency-domain  solvers  was 
signihcantly  reduced,  while  the  maximum  stable  CFL  of  the  split-domain  solver  re¬ 
mained  unchanged.  For  example,  for  the  a  =  0.3,  /  =  1.5  test  case  on  a  501  point 
grid,  the  split-domain  solver  was  able  to  compute  a  48-frequency  solution  with  a 
CFL  of  1.0,  the  maximum  CFL  for  the  MacCormack  discretization  scheme.  For  the 
same  conhguration,  the  direct,  time-domain,  and  frequency- domain  solvers  became 
unstable  when  the  CFL  exceeded  0.42. 

For  super-asymptotic  solutions,  increased  series  length  sometimes  had  a  detri¬ 
mental  effect  on  numerical  stability,  even  for  the  split-domain  solver.  All  sub- 
asymptotic  and  asymptotic  solutions,  however,  were  successfully  converged  to  a 
residual  (L2  norm)  of  5.0E-8  with  a  CFL  of  0.95,  independent  of  the  number  of 
frequencies  used. 

The  need  for  a  reduced  CFL  was  not  consistently  related  to  the  number  of 
super-asymptotic  frequencies  included  in  the  solution.  Only  test  cases  with  a  dis¬ 
turbance  frequency  of  3.0  required  a  reduced  CFL  to  obtain  48-frequency  solutions. 
In  these  cases,  reductions  of  40%  to  60%  in  CFL  were  required.  All  other  test  cases 
converged  with  an  un-reduced  CFL,  despite  the  fact  that  some  of  those  cases  had 
smaller  asymptotic  frequencies  than  some  of  the  /  =  3.0  cases. 
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3.4-3  Performance.  The  best  overall  run-time  performance  was  obtained 
from  the  split-domain  solver,  especially  for  large  N  solutions.  For  a  typical  49- 
frequency  solution,  the  split-domain  solver  ran  approximately  3  times  faster  (based 
on  CPU  time)  than  the  frequency-domain  and  time-domain  solvers,  and  7  times 
faster  than  the  direct  solver.  Factors  contributing  to  the  faster  run  time  were  the 
larger  time  step  allowed  by  the  higher  CFL,  and  the  fact  that  the  split-domain  scheme 
calculates  fluxes  in  the  time  domain  with  only  a  single  Fourier  transform  pair. 

A  third  factor  that  heavily  influenced  the  performance  of  all  the  multi-domain 
solvers,  including  the  split-domain  solver,  was  the  choice  of  iV  =  49.  As  discussed 
in  Section  2.2,  the  computational  cost  of  an  FFT,  and  thus  of  a  multi-domain  flux 
calculation,  is  highly  dependent  on  the  number  of  terms  being  transformed.  This  is 
dramatically  illustrated  by  examining  run  times  for  a  typical  48-frequency  solution. 

At  48  frequencies,  the  split-domain  solver  is  only  about  1.8  times  faster  than 
the  direct  solver,  while  the  time-domain  and  frequency  domain  solvers  are  about  2.4 
times  slower  than  the  direct  solver.  The  drop  in  performance  is  explained  by  looking 
at  the  number  of  terms  being  transformed  by  the  FFT  for  each  N.  For  N  =  49,  the 
number  of  terms  transformed  is  2N  -|-  1  =  99,  which  is  easily  factored  into  3^11^, 
and  thus  the  proportional  cost  of  the  FFT  is  close  to  the  best-case  N  log2  N.  But 
when  N  =  48,  the  number  of  terms  transformed  is  97,  which  is  prime  and  cannot 
be  factored  into  products  of  small  primes.  In  this  case,  the  cost  of  the  FFT  is 
proportional  to  the  worst-case  iV^.  As  a  result,  all  of  the  multi-domain  solvers  take 
appreciably  longer  to  run  when  N  is  reduced  from  49  to  48.  The  split-domain  solver 
manages  to  outperform  the  direct  solver,  but  only  because  it  is  running  at  a  higher 
CFL  and  computes  a  single  transform  pair  per  point  per  iteration. 

3. 5  Summary 

Large  amplitude,  time-periodic  solutions  to  Burgers’  equation  were  computed 
with  a  split-domain  harmonic  balance  solver,  and  compared  to  solutions  computed 
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using  prior  harmonic  balance  approaches.  The  split-domain  method  produced  so¬ 
lutions  comparable  to  those  of  the  prior  methods,  while  successfully  eliminating  a 
stability  restriction  experienced  by  those  methods  when  a  large  number  of  Fourier 
frequencies  are  included  in  the  solution.  This,  combined  with  a  reduction  in  the  num¬ 
ber  of  FFTs  required  to  implement  the  split-domain  method,  resulted  in  significantly 
reduced  run  times. 

Solutions  for  boundary  conditions  containing  moving  waves  ranging  from  smooth 
disturbances  to  strong  discontinuities  were  successfully  computed.  Comparison  with 
conventional  time-accurate  calculations  showed  that  the  error  in  the  harmonic  bal¬ 
ance  solutions  was  asymptotic  with  respect  to  the  number  of  frequencies  included  in 
the  approximating  solution.  When  the  number  of  frequencies  was  equal  to  or  greater 
than  the  asymptotic  frequency,  the  harmonic  balance  solutions  were  comparable  to 
the  time-accurate  solutions.  Several  factors  were  found  to  influence  the  asymptotic 
frequency,  including  disturbance  frequency,  the  strength  of  the  moving  wave,  and 
the  computational  grid  density. 
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IV.  Application  of  Adaptive  Split-Domain  Harmonic  Balance  to 

Euler’s  Equation 

4-1  Introduction 

This  chapter  documents  an  investigation  of  the  adaptive  harmonic  balance 
method  for  CFD.  The  primary  goal  of  the  investigation  was  to  establish  that  the 
energy-based  frequency  augmentation  approach  described  in  Section  2.5  reliably  and 
effectively  matches  harmonic  balance  frequency  content  to  local  flow  conditions, 
producing  accurate  solutions  in  less  time  than  a  non-adapted  harmonic  balance  ap¬ 
proach.  The  effects  of  grid  density,  augmentation  threshold,  pixelation  width,  and 
adaptation  scheduling  on  the  quality  of  adapted  solutions  and  performance  of  the 
adaptation  algorithm  were  examined.  Finally,  the  compatibility  of  the  adaptive  split- 
domain  harmonic  balance  approach  with  FAS  multigrid  convergence  acceleration  was 
conhrmed. 

To  accomplish  these  objectives,  adapted  harmonic  balance  solutions  were  com¬ 
puted  for  a  variety  of  supersonic  and  subsonic  inviscid  1-D  periodic  flows  governed  by 
the  1-D  Euler’s  equation.  Results  were  compared  with  non-adapted  harmonic  bal¬ 
ance  solutions,  and  with  solutions  obtained  using  conventional  time-accurate  CFD 
techniques. 


4-2  Solver  Implementation 

The  one-dimensional  Euler  equations  in  strong  conservation  form  are  given  by 


(9Q  5F 

dt  dx 


0 


(4.1) 
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where 


pu 

F  =  pu'^  +  p  (4.2) 

{Et  +  p)u 

and  p,  M,  and  p  are  the  density,  velocity,  and  pressure,  respectively.  Total  energy,  Et, 
is  dehned  as  p(e+  where  e  is  specihc  internal  energy.  All  quantities  are  nondi- 
mensional.  Together  with  the  perfect  gas  relation  and  an  assumption  of  a  constant 
ratio  of  specihc  heats,  7  =  1.4,  Eq.  4.1  represents  a  closed  system  of  equations. 

Equation  4.1  provides  the  form  of  the  time-domain  PDE  in  a  split-domain  har¬ 
monic  balance  implementation  of  the  1-D  Euler  equation.  In  the  harmonic  balance 
implementation,  the  dependent  vector  Q  is  replaced  by  a  new  vector  Q  which  is 
composed  of  the  2N  -|-  1  instances  of  the  vector  Q  sampled  at  times  f  =  (0,  At, 
2 At,  ■  ■  ■ ,  2N At),  with  At  =  ^(2W+i)  ’  period  of  oscillation  divided  by  the  number 
of  samples,  2N  -|-  1.  However,  because  the  equations  corresponding  to  each  sample 
are  independent  of  the  other  samples,  they  can  be  treated  independently.  For  this 
reason,  the  following  discussion  is  based  on  a  single  sample. 

A  cell-centered  hnite-volume  solver  (30)  was  written  to  solve  the  time-domain 
PDE.  The  hnite-volume  solver  is  based  on  a  discretization  of  the  integral  form  of 
Eq.  4.1 

^  f  QdV+  [  F-hdS  =  0.  (4.3) 

dr  Jv  Js 

In  Eq.  4.3,  dV  is  a  diherential  volume  element,  h  is  a  unit  vector  normal  to  the  control 
volume  surface,  and  dS  is  a  diherential  surface  element.  Integration  is  performed 
over  the  interior  (hrst  term)  and  surface  (second  term)  of  the  control  volume.  The 
physical  time  variable  t  has  been  replaced  by  a  pseudo-time  variable  r.  If  the  control 
volume  is  taken  to  be  a  single  grid  cell,  then  for  a  one-dimensional  computational 
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grid,  Eq.  4.3  reduces  to 


(4.4) 


+  -F'  =  0 

dr  *  * 

for  each  cell.  The  one-dimensional  cell  volume,  V,  is  just  the  cell  size  Ax,  and  Qi 
is  the  average,  or  cell  center,  value  of  the  conserved  variables  inside  the  cell.  The 
fluxes  at  the  left  and  right  faces,  F*  and  F*”,  are  constructed  by  averaging  the  fluxes 
evaluated  at  the  cell  centers  on  either  side  of  the  face.  For  a  grid  with  uniform  cell 
size,  this  scheme  is  second-order  accurate  in  space  and  is  equivalent  to  a  second-order 
central  difference  scheme. 

To  avoid  oscillations  at  discontinuities,  second  and  fourth  order  modihed  Jameson- 
Schmidt-Turkel  (JST)  artihcial  dissipation  (30,  31)  was  implemented.  JST  artihcial 
dissipation  adds  an  additional  term,  D,  of  the  form 

D  =  (D^  -  D^)Q  (4.5) 

to  the  left  hand  side  of  Eq.  4.4.  and  are  second  and  fourth  order  difference 
operators  dehned  as 


Qi  =  V(Ai+i/2  Q*  (4-6) 

=  (VA)(A,£f)VA)Q,  (4.7) 

where 

4+’i/2  =  A'(^)max(z/i,z/i+i)  (4.8) 

=  max(0,  77^2)  -  TsT^^Vi)  (4.9) 

Pi+i  -  2pi  +  pi_i 

Vi  =  -  .  (4.10) 

Pi+i  +  2pi  -h  Pi-1 


In  the  above  equations,  A  and  V  are  hrst-order  forward  and  backward  differ¬ 
ence  operators.  A*  =  |ui|  -|-  a*  is  the  maximum  eigenvalue  of  the  flux  Jacobian  matrix 
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1^,  Ai+1/2  =  l/2(Ai  +  Aj+i),  and  is  the  local  speed  of  sound.  The  parameters 
and  control  the  amount  of  dissipation  applied,  and  are  set  by  the  user  at  run 
time.  For  the  test  problems  described  below,  typical  values  of  ranged  from  0.4 
to  0.55,  while  values  of  ranged  from  0.004  to  0.04. 

Experience  with  the  code  showed  that  solution  quality  and  convergence  prop¬ 
erties  were  sometimes  improved  by  including  additional  second  order  dissipation. 
The  additional  dissipation  was  included  by  means  of  a  modihed  operator,  given 
by 

+  A,+i/2£® /2)AQ,.  (4.11) 

The  new  parameter  controls  the  amount  of  additional  second  order  dissipa¬ 
tion.  For  positive  second  order  dissipation  in  the  form  of  second  difference 

smoothing  is  applied  throughout  the  grid.  If  is  zero,  the  original  modihed  JST 
dissipation  is  recovered.  Typical  values  of  used  in  this  study  ranged  from  0.0 

for  supersonic  hows  to  0.001  for  subsonic  hows. 


Including  the  artihcial  dissipation  term  and  replacing  the  cell  volume  with  the 
cell  size.  Ax,  the  semi-discretized  one-dimensional  hnite-volume  Euler’s  equation 
becomes 


dQi 


D, 


(4.12) 


ch-  Ax'  • 

A  three-stage  Runge-Kutta  (RK)  scheme  (Eq.  4.13)  with  good  high-frequency  smooth¬ 
ing  properties  and  minimal  storage  requirements  (32)  was  used  to  advance  the  solu¬ 
tion  in  pseudo-time. 


Q*  =  Q”  +  0.35ArR(Q”) 

Q”  =  Q’^  +  0.6ArR(Q*)  (4.13) 

qn+i  _  Q-  +  ArR(Q**) 
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In  Eq.  4.13,  i?(Q)  is  the  right  hand  side  of  Eq.  4.12  evaluated  at  the  indicated  time 
levels.  Artihcial  dissipation  terms  were  evaluated  once  during  the  hrst  RK  stage, 
and  frozen  for  the  remaining  stages. 

Implementation  of  the  split-domain  harmonic  balance  approach  also  required 
the  solution  of  the  frequency-domain  ODE,  Eq.  2.23b.  Two  solution  approaches  were 
tested — an  exact  integration  approach,  given  by  Eq.  2.25,  and  a  numerical  approach 
based  on  the  Runge-Kutta  integration  scheme  applied  to  the  time-domain  PDE.  It 
was  found  that  the  exact  approach  provided  no  signihcant  improvement  in  solution 
quality  or  stability  compared  to  numerical  integration,  and  sometimes  slowed  or 
stalled  convergence  to  steady  state.  For  this  reason,  the  numerical  approach  was 
employed  for  the  results  documented  below. 

Boundary  conditions  were  enforced  by  setting  a  single  ghost  cell  outside  the 
boundaries  of  the  computational  domain  (see  Appendix  D).  Because  only  one  ghost 
cell  was  maintained,  the  fourth-order  dissipation  term,  Q,  was  set  to  zero  in 
the  cells  adjacent  to  the  boundaries.  All  boundary  conditions  were  calculated  and 
applied  in  the  time-domain. 

Local  time  stepping  was  employed  to  accelerate  convergence  to  steady  state. 
The  maximum  stable  time  step  allowed  in  the  cell  was  determined  as  follows: 

^Tmax  =  CEL  min  •  (4.14) 

samples  \  / 

Local  time  steps  were  calculated  every  iteration  in  conjunction  with  the  artihcial 
dissipation  term,  and  used  immediately.  This  meant  that  the  time  step  used  for 
the  hrst  frequency-domain  ODE  integration  (step  one  of  the  split-domain  iteration. 
Section  2.4)  was  diherent  from  that  used  for  steps  three  and  hve,  the  time-domain  and 
second  frequency-domain  integration.  This  did  not  ahect  the  steady-state  solution, 
however,  because  the  change  in  time  step  from  one  iteration  to  another  became 
negligible  as  the  solution  approached  steady  state. 
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The  maximum  stable  CFL  for  the  solver  was  1.7.  In  a  number  of  test  cases  with 
large  sample  rates,  a  reduced  CFL  was  required.  A  simple  variable  CFL  function 
was  implemented  that  linearly  varied  the  CFL  from  a  user-specihed  maximum  at  the 
lowest  sample  rate  to  a  user-specihed  minimum  at  the  highest  sample  rate.  During 
the  early  stages  of  the  solution  when  sample  rates  were  low,  a  large  CFL  was  applied. 
As  the  solution  progressed  and  frequencies  were  added  to  the  solution,  the  CFL  was 
reduced  to  maintain  a  stability  solution.  The  linear  scaling  approach  provided  a 
conservative  but  robust  scaling  of  CFL  with  sample  rate. 

In  addition  to  local  time  stepping,  full  approximation  storage  (FAS)  multigrid 
convergence  acceleration  (28,  33)  was  implemented  as  described  in  Section  2.6.  Two 
different  restriction  operators  were  used  to  transfer  the  hne-grid  solution  to  the 
next  coarsest  grid  (33).  Split-domain  residuals  (Eq.  2.39)  were  transfered  using  a 
conservative  transfer  operator  dehned  by 

ILiRm  =  T  5^(n+ifls+i).  (4.15) 

Vk 

where  the  snbscripts  k  and  k  +  1  identify  coarse  and  hne  grid  valnes  respectively, 
and  V  is  cell  volume.  Summation  occnrs  over  the  two  hne  grid  cells  making  np  each 
coarse  grid  cell.  A  volume-weighted  transfer  operator,  dehned  by 

rfc  A  _  ^(^fc+lQfc+l)  //I  1 

was  used  to  transfer  the  solution  valnes  themselves.  Immediately  after  restriction, 
initial  local  time  steps  for  the  coarse  grid  were  calculated  according  to  Eq.  4.14. 
Once  a  coarse  grid  correction  was  computed,  it  was  transfered  to  the  next  hnest  grid 
nsing  linear  interpolation  as  shown  in  Fig.  4.1. 
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Qfe+ibi  —  0.75  Qfe|i  +  0.25  Qfe|i+i 
Qfc+ibi+i  =  0.25  Qfc|*  + 0.75  Qfc  |i+l 

Figure  4.1  Linear  Interpolation  from  Coarse  Grid  to  Fine  Grid 
4-3  Test  Configurations 

Adaptive  frequency  augmentation  was  applied  to  several  one-dimensional  test 
cases.  The  first  set  of  cases  consisted  of  a  family  of  flow  fields  with  sinusoidally 
varying  supersonic  inflow.  For  these  cases,  the  inflow  Mach  number  was  varied 
according  to 

M  =  2.0 -|- a  sin  (cut)  (4-17) 

while  maintaining  a  constant  nondimensional  density  and  static  pressure  of  1.0  and 
0.17857,  respectively.  The  disturbance  frequency,  oj  was  chosen  so  that  approxi¬ 
mately  two  complete  disturbance  cycles  would  occur  each  flow-through  period  on  a 
unit  grid.  The  magnitude  of  the  sinusoidal  disturbance  was  varied  to  achieve  different 
flow  characteristics. 

Pressure  distributions  at  several  snapshots  in  time  for  a  flow  held  with  a  =  0.25 
are  shown  in  Fig.  4.2.  These  solutions  were  obtained  from  a  fully-developed  time- 
accurate  Roe  solver  calculation.  Because  the  supersonic  how  was  hyperbolic  in  time 
with  all  how  information  traveling  in  the  downstream  direction,  the  fully-developed 
solution  was  achieved  very  quickly.  The  snapshots  plotted  in  Fig.  4.2  were  take  after 
5  1/2  how-through  periods. 
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Figure  4.3  Flow  Interactions  for  Simulated  Oscillating  Piston  in  an  Open  Tube 

The  last  test  case  consisted  of  a  purely  subsonic  flow  field  loosely  approximating 
the  flow  in  an  open  tube  with  an  oscillating  piston  at  one  end.  The  piston  action 
was  simulated  by  imposing  a  sinusoidal  velocity  u(t)  =  asin(ct;t)  at  a  fixed  boundary 
while  maintaining  a  zero  pressure  gradient  and  constant  total  enthalpy.  The  open  end 
of  the  tube  was  modeled  with  a  characteristic-variable  far-held  boundary  condition 
(34).  Details  of  the  boundary  condition  implementation  are  given  in  Appendix  D. 

While  these  boundary  conditions  do  not  exactly  model  a  physical  system,  they 
do  result  in  a  complex,  periodic,  subsonic  how  held  containing  alternating  right¬ 
running  shock  and  expansion  waves  interacting  with  rehections  from  the  open  end 
of  the  tube  (Fig.  4.3).  Initial  how  conditions  of  u  =  0,  p  =  1,  and  p  =  0.7142857 
were  assumed.  By  varying  the  length  of  the  tube  and  the  magnitude  and  frequency 
of  the  imposed  sinusoidal  piston  velocity,  the  how  held  was  tuned  so  that  a  sta¬ 
tionary  periodic  how  held  with  moderately  strong  features  was  produced.  Tuning 
was  accomplished  by  computing  a  large  number  of  piston  cycles  with  a  validated 
time-accurate  CFD  solver,  starting  from  a  zero-velocity  initial  solution.  The  time 
history  of  the  mass  in  the  tube  was  used  to  determine  when  a  stationary  solution  was 
achieved.  The  chosen  conhguration  achieved  fully  developed  how  in  approximately 
150  piston  cycles,  as  shown  in  Fig.  4.4. 

Figure  4.5  contains  pressure  plots  for  the  fully  developed  subsonic  test  how 
held  at  several  points  in  time,  computed  during  the  199th  piston  cycle  of  the  time- 
accurate  calculation. 
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Figure  4.4  Time-accurate  Development  of  the  Stationary  Flow  Field  for  the 
Piston-In-Tube  Configuration 

All  harmonic  balance  solutions  for  both  the  supersonic  and  subsonic  test  cases 
were  converged  to  an  overall  residual  (Eq.  2.29),  of  l.Oe-9,  or  approximately  6.5  or¬ 
ders  of  magnitude.  For  comparison  purposes,  time-accurate  solutions  such  as  those 
discussed  above  were  also  computed.  A  summary  of  case  names  and  dehning  param¬ 
eters  of  all  the  test  cases  is  contained  in  Table  4.1. 

Table  4.1  Summary  of  Test  Conhguration  Parameters 
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4-4  Tuning  the  Adaptive  Solver 

The  following  paragraphs  discuss  several  user-dehnable  parameters  that  con¬ 
trol  frequency  augmentation  and  their  impact  on  solution  accuracy  and  solver  per¬ 
formance.  For  a  detailed  discussion  of  the  adaptive  algorithm  and  these  parameters, 
see  Section  2.5. 

4-4-i  Selection  of  Augmentation  Threshold.  The  key  parameter  of  the 
adaptive  harmonic  balance  method  is  the  augmentation  threshold,  E thresh-  This 
threshold  value  indirectly  controls  the  sample  rate  distribution,  and  thus  the  run 
time  of  a  calculation  and  the  accuracy  of  the  resulting  solution. 

To  study  the  impact  of  Ethresh  on  frequency  content,  accuracy,  and  run  time, 
solutions  were  generated  for  both  supersonic  and  subsonic  conhgurations  with  thresh¬ 
olds  ranging  from  5.0e-4  to  l.Oe-10.  Examples  of  the  frequency  content  and  solution 
quality  resulting  from  different  augmentation  thresholds  are  shown  in  Figs.  4. 6-4. 9. 
These  hgures  show  frequency  and  pressure  distributions  for  conhgurations  SS3  and 
SUl,  respectively,  for  Ethresh  =  l.Oe-4,  l.Oe-7,  and  l.Oe-10.  The  pressure  distributions 
were  reconstructed  from  the  computed  Fourier  coefficients  at  time  t  =  0  relative  to 
the  disturbance  period.  It  is  important  to  note  that  the  reconstruction  could  just  as 
easily  have  been  generated  for  any  time  with  equal  hdelity. 

Figure  4.6  shows  adapted  frequency  distributions  for  conhguration  SS3.  From 
this  plot,  one  can  determine  the  number  of  Fourier  frequencies  included  at  each  cell  in 
the  computational  grid  for  each  of  the  three  augmentation  thresholds.  As  seen  in  the 
hgure,  there  was  a  signihcant  increase  and  rehnement  of  the  frequency  distribution 
with  the  change  from  l.Oe-4  to  l.Oe-7.  The  increased  frequency  content  was  accompa¬ 
nied  by  a  signihcant  improvement  in  solution  quality  as  shown  in  Fig.  4.7.  However, 
further  decreasing  the  threshold  to  l.Oe-10  resulted  in  no  perceptible  change  to  the 
reconstructed  pressures,  despite  another  large  increase  in  frequency  content. 


4-12 


Figure  4.6  Adapted  Frequency  Distributions  for  Three  Different  Augmentation 
Thresholds,  Case  SS3,  500  cells 


Figure  4.7  Nondimensional  Static  Pressure  for  Adapted  Solutions  at  Three  Differ 
ent  Augmentation  Thresholds,  Case  SS3,  500  Cells 


The  adapted  frequency  content  of  the  subsonic  case  also  increased  significantly 
with  each  decrease  in  Ethresh  (Fig.  4.8).  Unlike  the  supersonic  case,  however,  the 
subsonic  pressure  distributions  for  thresholds  of  l.Oe-4  and  l.Oe-7  are  nearly  identical 
(Fig.  4.9).  The  increased  freqnency  content  resnlting  from  the  decrease  in  threshold 
did  not  improve  the  solution.  A  further  threshold  decrease  to  l.Oe-10  resulted  in 
a  reduction  in  solution  quality,  as  high-frequency  oscillations  appeared  upstream  of 
the  strongest  shocks. 

The  resnlts  for  both  the  snbsonic  and  snpersonic  cases  are  consistent  with  the 
Burgers’  equation  results  presented  in  Section  3.4.  There  it  was  found  that  including 
more  frequencies  in  the  harmonic  balance  solution  did  not  always  improve  the  so- 
Intion;  there  was  an  asymptotic  freqnency  beyond  which  no  improvement  occurred. 
The  same  should  be  true  for  augmentation  threshold,  since  it  indirectly  controls 
frequency  content.  For  the  current  test  conhgurations,  solution  qnality  did  not  im¬ 
prove  as  Ethresh  was  decreased  below  5.0e-8  (supersonic)  and  l.Oe-4  (subsonic).  In 
Section  3.4.2,  it  was  also  found  that  in  some  cases,  a  reduced  CFL  was  required  to 
obtain  a  stable  solution  incorporating  frequencies  above  the  asymptotic  frequency. 
This  was  the  case  for  the  snbsonic  conhgnration  with  thresholds  below  l.Oe-6  (Ta¬ 
ble  4.2).  Evidence  suggests  that  the  low  CFL  required  to  obtain  a  solution  with 
Ethresh  =  l.Oe-10  contributed  to  the  appearance  of  high-freqnency  oscillations  in  the 
subsonic  solution  (see  Appendix  B). 

Table  4.2  Maximum  Fourier  Frequency  and  Stable  CFL  for  Configuration  SUl, 
500  Cells,  for  Decreases  in  Augmentation  Threshold 


Threshold 

Max  N 

Max  CFL 

l.Oe-6 

40 

1.7 

l.Oe-7 

45 

1.5 

l.Oe-8 

58 

1.1 

l.Oe-9 

67“ 

1.0 

l.Oe-10 

82“ 

0.7 

“  Maximum  allowed  for  run 
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The  effect  of  Ethresh  on  problem  size  (given  by  the  average  frequency  content 
of  the  solution),  and  run  time  is  shown  in  Figs.  4.10  and  4.11.  Figure  4.10  shows 
the  average  number  of  Fourier  frequencies  included  in  each  adapted  solution.  For 
both  the  supersonic  and  subsonic  conhgurations,  the  average  frequency  content  was 
inversely  proportional  to  the  logiQ^E thresh)-  The  size  of  the  supersonic  adapted  har¬ 
monic  balance  problem  grew  approximately  8  times  from  the  largest  threshold  to  the 
smallest,  while  the  size  of  the  subsonic  problem  grew  approximately  3.5  times. 

The  growth  in  run  time  with  decreasing  Ethresh  is  shown  in  Fig.  4.11.  Run  time 
for  the  supersonic  case  grew  at  approximately  the  same  rate  as  the  problem  size, 
increasing  6.7  times  from  the  largest  to  smallest  thresholds.  Down  to  a  threshold 
of  5e-8,  run  time  for  the  subsonic  case  also  grew  at  the  same  rate  as  problem  size. 
For  smaller  thresholds,  however,  the  subsonic  run  time  grew  at  approximately  twice 
the  rate  of  problem  size.  This  was  due  to  the  reduced  CFL  required  at  these  low 
thresholds. 


Figure  4.10  Change  in  Average  Frequency  Content  with  Changing  Augmentation 
Threshold  for  Conhgurations  SS3  and  SUl  on  500  Cell  Grid 
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Figure  4.11  Change  in  Relative  Compute  Time  with  Changing  Augmentation 
Threshold  for  Conhgurations  SS3  and  SUl  on  500  Cell  Grid 

These  results  suggest  that  a  reasonable  lower  limit  for  the  augmentation  thresh¬ 
old  would  be  5.0e-8.  Adapting  to  a  smaller  threshold  increased  run  times,  but  re¬ 
sulted  in  no  improvement  in  solution  quality.  Determining  a  suitable  upper  limit  for 
the  augmentation  threshold  is  less  clear.  Based  on  these  results,  a  conservative  upper 
bound  for  a  high-fidelity  solution  would  be  on  the  order  of  5.0e-7.  A  larger  thresh¬ 
old  may  be  acceptable  for  some  problems,  however,  as  was  seen  with  configuration 
SUl,  where  a  threshold  of  l.Oe-4  produced  comparable  results  (c.f..  Fig.  4.9).  It  is 
likely  that  a  maximum  acceptable  value  will  vary  widely  depending  on  the  flow  being 
modeled,  the  grid  used,  and  the  desired  solution  quality.  If  a  lower-fidelity  solution 
is  acceptable  (e.g.,  only  the  time-average  solution  is  of  interest)  then  a  much  higher 
threshold  would  suffice. 

4 -4 Fringe  Width.  None  of  the  test  cases  required  fringe  augmentation. 
Therefore  the  fringe  augmentation  width  was  set  to  zero. 
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4-4-3  Pixelation  Width.  Solutions  for  both  subsonic  and  supersonic  test 
cases  were  insensitive  to  the  adaptation  pixelation  width.  Widths  as  small  as  1 
(no  pixelation)  and  as  large  as  25  were  applied,  with  no  signihcant  change  in  the 
solntion.  For  the  snbsonic  case,  solutions  with  larger  pixelation  widths  converged 
slightly  faster  than  those  with  smaller  widths  (4%  faster  for  a  width  of  25,  vice  a 
width  of  1);  but  for  the  supersonic  case  no  consistent  advantage  was  observed. 

Based  on  these  resnlts,  there  is  no  reqnirement  for  freqnency  pixelation.  A 
slight  performance  improvement  might  be  realized  by  inclnding  pixelation,  but  any 
advantage  will  be  case  dependent.  Other  factors,  such  as  multigrid  implementation, 
may  increase  the  importance  of  pixelation,  however  (see  Section  2.6). 

4-4-4  Adaptation  Trigger.  Solutions  were  computed  for  various  values  of 
the  primary  residual-based  adaptation  triggers,  labeled  ki  (initial)  and  K2  (subse¬ 
quent).  To  ensnre  that  augmentation  was  based  on  these  triggers,  the  iteration-based 
triggers  were  disabled. 

The  adaptation  triggers  were  fonnd  to  have  little  impact  on  the  hnal  solution. 
They  did,  however,  have  a  signihcant  impact  on  run  times.  The  behavior  for  a  given 
set  of  triggers  was  qnite  different  for  supersonic  and  subsonic  solutions.  For  the 
supersonic  conhgnrations,  tests  showed  that  it  was  best  to  adapt  after  only  a  small 
reduction  in  residual.  The  hyperbolic  nature  of  the  supersonic  how  helds  meant  that 
any  disturbance  introduced  in  the  how  had  to  propagate  ont  of  the  grid  before  the 
solution  converged.  Thus,  there  was  no  advantage  in  converging  a  solution  prior  to 
adding  new  frequencies.  The  optimum  trigger  values  for  the  supersonic  cases  were 
found  to  be  Ki  =  0.25  and  K2  =  0.  In  comparison,  a  solution  with  =  2.0  and 
fi;2  =  0.1  took  50%  longer  to  compute. 

The  opposite  behavior  was  observed  for  the  subsonic  case.  Here,  the  elliptic 
nature  of  the  problem  meant  that  converging  the  low-frequency  solutions  before 
adapting  provided  a  distinct  performance  advantage.  The  best  run  time  for  the 
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subsonic  case  were  achieved  with  ki  =  2.5  and  K2  =  0.1.  The  longest  run  time 
(approximately  55%  longer)  was  achieved  with  ki  =  0.5  and  K2  =  0.1.  Attempts  to 
run  with  ki  =  0.25  resulted  in  unstable  solutions. 

One  subsonic  test,  conducted  with  ki=  2.75,  never  adapted.  The  solution 
residual  at  the  initial  frequency  content  {N  =  3)  plateaued  after  a  drop  of  approxi¬ 
mately  2.55  orders  of  magnitude,  and  failed  to  achieve  the  specihed  drop.  Had  the 
iteration-based  triggers  been  enabled,  they  would  have  ensured  that  the  solution  was 
eventually  adapted. 

4.5  Results 

4-5.1  Accuracy.  Figures  4.12-4.15  compare  results  for  both  adapted  and 
non-adapted  harmonic  balance  solutions  and  a  conventional  time-accurate  solution. 
Figs.  4.12  and  4.13  contain  pressure  and  coefficient  magnitudes  for  conhguration  SS4, 
while  Figs.  4.14  and  4.15  contain  similar  information  for  configuration  SUl. 

The  pressure  distributions  shown  for  conhguration  SS4  (Fig.  4.12)  are  typical 
of  all  the  supersonic  conhgurations.  The  reconstructed  pressures  for  the  adapted  and 
non-adapted  solutions  are  equivalent,  and  compare  favorably  with  the  time-accurate 
calculation.  Both  include  slight  second  order  Gibbs  effects  at  the  shocks  that  would 
be  removed  with  an  increase  in  artihcial  dissipation.  Examining  the  magnitudes 
of  the  momentum  term  Fourier  coefficients  (Fig.  4.13)  conhrms  that,  except  for 
some  small  differences  in  the  highest  frequencies  and  near  frequency  transitions, 
the  adapted  and  non-adapted  results  are  essentially  identical.  This  suggests  that 
the  omitted  frequencies  in  the  augmented  solution  have  no  signihcant  impact  on 
accuracy. 

The  reconstructed  pressure  distributions  for  the  adapted  and  non-adapted  har¬ 
monic  balance  solutions  for  conhguration  SUl  (Figs.  4.14)  are  comparable.  There  are 
slight  differences  in  the  two  solutions,  however,  as  can  be  seen  by  examining  the  co¬ 
efficient  magnitude  contours  in  Fig.  4.15.  Unlike  the  supersonic  case,  there  are  small 
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Figure  4.12  Computed  Nondimensional  Static  Pressure  for  Adapted  (Augmenta¬ 
tion  Threshold  7.0e-8)  and  Non-adapted,  49-Frequency  Harmonic  Bal¬ 
ance  Solutions,  Configuration  SS4,  601  Cells,  Compared  with  1000 
Cell  Time-accurate  Solution 


Figure  4.13  Fourier  Coefficient  Magnitudes  (Momentum  Component)  for  Adapted 
(Augmentation  Threshold  7.0e-8)  and  Non-adapted,  49-Frequency 
Harmonic  Balance  Solutions,  Configuration  SS4,  601  Cells 
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Figure  4.14  Computed  Nondimensional  Static  Pressure  for  Adapted  (Augmenta¬ 
tion  Threshold  7.0e-8,  Variable  CFL)  and  Non-adapted,  58-Frequency 
Harmonic  Balance  Solutions,  Conhguration  SUl,  601  Cells,  Compared 
with  1000  Cell  Time-accurate  Solution 


Figure  4.15  Fourier  Coefficient  Magnitudes  (Momentum  Component)  for  Adapted 
(Augmentation  Threshold  7.0e-8,  Variable  CFL)  and  Non-adapted, 
58-Frequency  Harmonic  Balance  Solutions,  Conhguration  SUl,  601 
Grid  Cells 


4-21 


but  noticeable  differences  between  the  adapted  and  non-adapted  harmonic  balance 
Fourier  coefficients.  These  differences  result  from  the  use  of  a  variable  CFL  in  the 
adapted  solution.  The  steady-state  solution  computed  by  the  split-domain  harmonic 
balance  method  has  a  small  dependence  on  time  step  size  (see  Appendix  B).  Appli¬ 
cation  of  CFL  scaling  resulted  in  larger  time  steps  in  some  cells,  and  thus  slightly 
changed  the  steady-state  solution.  The  two  solutions  would  be  made  nearly  iden¬ 
tical  by  freezing  the  CFL  of  the  adapted  solution  at  the  same  value  used  for  the 
non-adapted  solution,  but  this  would  increase  run  time  and  result  in  only  minute 
changes  to  the  reconstructed  pressures. 

Both  adapted  and  non-adapted  harmonic  balance  solutions  for  conhguration 
SUl  show  some  discrepancies  when  compared  to  the  time-accurate  solution,  par¬ 
ticularly  in  the  trough  regions.  These  discrepancies  were  traced  to  the  boundary 
condition  at  the  open  (right)  end  of  the  grid.  A  comparison  of  the  computed  pres¬ 
sure  at  that  boundary  for  the  harmonic  balance  and  time  accurate  solutions  is  shown 
in  Fig.  4.16.  The  pressure  recovery  immediately  following  the  period  of  negative  ve¬ 
locity  (inflow)  is  much  sharper  for  the  time-accurate  solution  than  for  the  harmonic 
balance  solution.  The  steady-state  solutions  calculated  for  time  samples  in  the  re¬ 
covery  region  have  a  shock  wave  located  at  the  exit  boundary  (Fig.  4.17).  The  far 
held  boundary  condition  applied  at  that  boundary  was  not  designed  for  such  an 
extreme  gradient;  the  pressures  at  the  boundary  were  over-predicted.  The  same 
boundary  condition  was  used  in  the  time-accurate  code,  but  in  that  case,  the  shock 
at  the  boundary  was  moving,  and  so  its  effect  on  the  boundary  condition  was  less 
severe.  The  end  result  was  a  delay  and  attenuation  of  the  left-running  expansion 
wave  rehected  from  the  right  boundary. 

4-5.2  Grid  Density.  The  maximum  signihcant  Fourier  frequency  for  any 
harmonic  balance  solution  is  highly  dependent  on  grid  density.  High  frequency  co¬ 
efficients  vary  rapidly  in  the  spatial  dimensions,  and  require  a  hne  grid  for  accurate 
resolution.  For  a  grid  too  coarse  to  support  a  given  Fourier  frequency,  the  coeffi- 
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Figure  4.16  Computed  Pressure  for  Time-accurate  and  Harmonic  Balance  Solu¬ 
tions  at  the  Open  End  of  the  Tube 


Figure  4.17  Steady-state  Harmonic  Balance  Pressure  Distribution  at  Time  Sample 
Just  after  Period  of  Negative  Velocity 

cient  magnitudes  are  damped,  and  less  energy  is  contained  in  that  frequency.  This 
establishes  a  natural  ceiling  for  frequency  augmentation,  which  leads  to  an  observed 
beneht  of  the  adaptive  harmonic  balance  approach  -  it  automatically  matches  the 
included  frequency  content  to  the  computational  grid.  The  variation  in  adapted 
frequency  distribution  with  grid  density  is  illustrated  in  Figs.  4.18  and  4.19,  which 
show  the  hnal  frequency  distributions  at  several  grid  densities  for  conhgurations  SS3 
and  SUl,  respectively. 
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Figure  4.18  Variation  in  Adapted  Frequency  Distributions  with  Changing  Grid 
Density  for  Test  Conhguration  SS3,  Augmentation  Threshold  5.0e-8 


Figure  4.19  Variation  in  Adapted  Frequency  Distributions  With  Changing  Grid 
Density  for  Test  Conhguration  SUl,  Augmentation  Threshold  5.0e-8 

4.5.3  Performance.  To  gauge  the  performance  beneht  of  the  adaptive 
harmonic  balance  approach,  two  performance  metrics  were  examined.  The  hrst 
metric  measured  problem  size  reduction,  as  measured  by  the  difference  between  the 
maximum  and  average  frequency  content  of  an  adapted  solution.  The  second  metric 
was  the  run  time  required  to  converge  a  solution  to  steady  state. 
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For  each  of  the  test  conhgurations  SS2,  SS3,  SS4,  and  SUl,  an  adapted  solution 
was  calculated  with  an  augmentation  threshold  of  7.0e-8,  optimal  triggers,  and  a 
pixelation  width  of  7.  No  multigrid  acceleration  was  used.  A  non-adapted  solution 
was  then  calculated,  again  with  no  multigrid,  based  on  the  highest  frequency  content 
in  the  adapted  solution.  Both  adapted  and  non-adapted  solutions  for  conhgurations 
SS2,  SS3,  and  SS4  were  computed  with  a  hxed  CFL  of  1.7.  For  case  SUl,  the  adapted 
solution  was  computed  with  a  CFL  that  scaled  from  1.7  at  the  lowest  frequency  to 
1.35  at  the  highest  frequency,  while  the  non-adapted  solution  required  a  hxed  CFL 
of  1.3. 

The  reduction  of  frequency  content  and  run  time  for  the  adapted  solution 
is  shown  in  Fig.  4.20.  In  all  cases,  adaptation  resulted  in  reduced  run  times,  de¬ 
spite  relatively  small  reductions  in  frequency  content.  The  subsonic  case  showed  the 
most  improvement,  with  more  than  a  50%  reduction  in  run  time  compared  to  the 
non-adapted  case.  Much  of  this  reduction  was  due  to  the  variable  CFL.  The  higher 
average  CFL,  coupled  with  the  efficiency  gained  by  solving  the  large-scale  how  struc¬ 
tures  with  fewer  frequencies,  resulted  in  a  run  time  reduction  20%  larger  than  the 
problem  size  reduction. 

Since  the  adapted  and  non-adapted  supersonic  cases  were  computed  with  the 
same  CFL,  the  adapted  solutions  had  no  time-step  advantage.  Due  to  the  overhead 
of  the  adaptive  approach,  the  number  of  iterations  required  to  converge  the  adapted 
solutions  increased  relative  to  the  non-adapted  solutions  by  as  much  as  15%.  Despite 
this  increase  in  iterations,  the  reduced  average  frequency  content  still  resulted  in  a 
small  but  signihcant  reduction  in  run  time. 

4-6  Multigrid 

To  demonstrate  that  the  adaptive  harmonic  balance  technique  is  compatible 
with  multigrid  acceleration,  FAS  multigrid  was  applied  to  one  test  case.  A  low- 
disturbance-magnitude  supersonic  case  (conhguration  SSI)  was  solved  on  a  1025 
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Figure  4.20  Reduction  in  Average  Frequency  Content  and  Run  Time  for  Adapted 
vs  Non-adapted  Harmonic  Balance  Solutions 

point  grid  using  a  5-level  FAS  Multigrid  V  cycle.  This  configuration  resulted  in  a 
solution  that  is  relatively  smooth  over  most  of  the  grid.  The  adapted  frequency 
distribution  contained  just  2  frequencies  at  the  inlet,  and  gradually  increased  in  1- 
frequency  increments  to  7  frequencies  at  the  exit.  A  representative  pressure  plot  for 
this  conhguration  computed  with  and  without  multigrid  convergence  acceleration  is 
shown  in  Fig.  4.21. 

One  minor  modihcation  to  the  augmentation  algorithm  was  made  to  simplify 
implementation  of  the  multigrid  scheme.  The  primary  residual-based  adaptation 
trigger  was  replaced  with  a  trigger  based  on  the  number  of  multigrid  cycles  com¬ 
pleted.  For  this  conhguration,  adapting  every  3  multigrid  cycles  was  found  to  be 
effective. 

Run  times  for  both  adapted  and  non-adapted  solutions,  with  and  without 
multigrid,  are  given  in  Table  4.3.  In  both  the  adapted  and  non-adapted  cases,  run 
times  with  multigrid  acceleration  are  approximately  half  those  without.  The  non- 
adapted  solution  benehts  slightly  more  than  the  adapted  solution. 
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Figure  4.21  Nondimensional  Static  Pressure  for  Adapted  (Augmentation  Thresh¬ 
old  7.0e-8)  Harmonic  Balance  Solutions,  Conhguration  SSI,  1025 
Cells,  With  and  Without  Multigrid  Acceleration 


Table  4.3  Run-time  Performance  Comparison  for  Adapted  vs  Non-adapted  Har¬ 
monic  Balance  Solutions,  With  and  Without  Multigrid  Acceleration, 
Case  SSI,  1025  Cells,  in  Seconds 


Non-Multigrid 

Multigrid 

7  Freq 

Adapted 

7  Freq 

Adapted 

102.3 

88.4 

50.9 

45.7 

A  7  Summary 

The  adaptive  split-domain  harmonic  balance  method  was  successfully  applied 
to  a  variety  of  supersonic  and  subsonic  one-dimensional  flow  helds  containing  strong 
moving  shocks.  The  energy-based  augmentation  approach  reliably  identihed  cells 
where  additional  frequency  content  was  needed  and  could  be  supported  by  the  com¬ 
putational  grid.  The  resulting  adapted  harmonic  balance  solutions  were  equivalent 
to  non-adapted  harmonic  balance  solutions,  and  compared  well  with  conventional 
time-accurate  solutions. 
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Adaptation  scheduling  was  found  to  have  a  signihcant  impact  on  the  run-time 
performance  of  the  adaptive  solver.  In  some  cases,  applying  the  wrong  scheduling 
approach  more  than  doubled  solution  compute  time.  For  supersonic  flows,  rapid 
adaptation  with  minimal  flow  development  time  produced  the  best  performance.  For 
subsonic  flows,  a  scheduling  approach  that  allowed  more  time  for  flow  development 
between  adaptations  was  best. 

The  compatibility  of  the  adaptive  split-domain  harmonic  balance  approach  and 
FAS  multigrid  acceleration  was  demonstrated.  Solutions  were  computed  with  and 
without  frequency  augmentation,  and  with  and  without  multigrid  acceleration.  The 
relative  performance  benefits  of  each  approach  remained  consistent. 

For  all  test  conhgurations,  adapted  harmonic  balance  solutions  took  less  time 
to  compute  than  equivalent  non-adapted  solutions.  Reductions  from  25%  to  greater 
than  50%  were  observed.  These  reductions  were  obtained  despite  the  fact  that  most 
of  the  test  configurations  contained  strong  discontinuities  throughout  the  solution 
domain,  and  thus  had  high  average  adapted  frequency  content. 
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V.  Application  of  Adaptive  Split-Domain  Harmonic  Balance  to  an 

Unsteady  Diverging  Nozzle 

5. 1  Introduction 

The  test  problems  solved  in  the  previous  chapter  using  the  frequency  augmen¬ 
tation  approach  were  useful  for  developing  the  adaptive  harmonic  balance  approach, 
but  they  were  not  ideal  for  demonstrating  performance  improvements.  In  those  test 
problems,  there  were  moving  shocks  throughout  the  computational  domain,  which 
led  to  a  high  average  frequency  content.  In  a  transonic  turbomachinery  application, 
one  would  expect  large  variations  in  flow  time  response  at  different  points  in  the 
computational  domain,  from  near  steady-state  flow  far  upstream  of  a  blade  row,  to 
a  region  of  strong  moving  shocks  near  the  leading  edge  of  the  blades. 

The  purpose  of  the  analysis  documented  in  this  chapter  is  to  demonstrate  the 
adaptive  split-domain  harmonic  balance  approach  on  a  problem  representative  of 
a  transonic  turbomachinery  problem.  The  chosen  test  case  consists  of  a  diverging 
nozzle  with  supersonic  inflow  and  subsonic  outflow.  Unsteadiness  is  introduced  by 
varying  the  subsonic  outflow  conditions.  The  resulting  flow  held  contains  steady- 
state  how  over  most  of  the  hrst  half  of  the  nozzle,  a  narrow  region  near  the  center 
of  the  nozzle  with  an  oscillating  normal  shock,  and  smoothly  varying  how  over  the 
latter  part  of  the  nozzle.  These  how  features  make  a  good  test  case  for  the  frequency 
augmentation  algorithm.  The  case  contains  mostly  smooth  how  with  low  (or  no) 
frequency  content.  It  abruptly  transitions  to  a  small  region  where  the  time  response 
of  the  how  is  essentially  a  square  wave,  and  thus  has  high  frequency  content.  It  then 
abruptly  transitions  back  to  low  frequency  content.  The  robustness  of  the  algorithm 
is  demonstrated  by  handling  the  transitions,  and  the  ehectiveness  is  demonstrated 
by  a  low  average  frequency  content  and  reduced  computation  time. 
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5.2  Solver  Implementation 

The  one- dimensional  split-domain  harmonic  balance  Euler  solver  developed  in 
Chapter  IV  was  modihed  to  solve  the  quasi-l-D  Euler  equation.  The  quasi- 1-D  Euler 
equation  (35)  in  hnite  volume  form  is  given  by 
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where  A  is  the  cross  sectional  area  of  the  nozzle  and  the  other  quantities  are  as 
dehned  for  the  1-D  Euler  equation  in  Section  4.2.  Substituting  the  approximation 
^  ^  and  assuming  the  computational  grid  does  not  change  with  time,  Eq.  5.1 

reduces  to 

+  A'F;  -  A\¥\  -  p,{A\  -  y )  =  0  (5.5) 

for  the  ith  grid  cell,  where  A  is  the  average  area  over  the  cell,  and  A^  and  A^  are 
the  cross  sectional  areas  at  the  left  and  right  cell  interfaces  respectively.  Adding 


5-2 


artificial  dissipation  as  discussed  in  Section  4.2,  and  recognizing  that  AiAxi  =  Vi, 
the  approximate  volume  of  a  vertical  slice  through  the  nozzle,  the  semi-discretized 
quasi- 1-D  Euler  equation  becomes 

^  -  F‘4  -  p,(A'  -  4)  +  D.).  (5.6) 

Modifying  the  1-D  Euler  code  to  solve  the  Quasi-l-D  Euler  equations  thus  required 
changing  the  dehnition  of  the  cell  volume,  multiplication  of  the  left  and  right  fluxes 
by  the  appropriate  cross  sectional  areas,  and  inclusion  of  the  pressure  term.  A 
characteristic-based,  specihed-density  outflow  boundary  condition  (35)  designed  for 
the  diverging  nozzle  was  also  incorporated  (See  Appendix  D).  Other  aspects  of  the 
solver  remain  as  described  in  Section  4.2. 


5.3  Test  Configuration 


The  quasi- 1-D  adaptive  split-domain  harmonic  balance  solver  was  applied  to 
the  problem  of  flow  through  a  diverging  nozzle  with  constant  supersonic  inflow  and 
unsteady  subsonic  outflow  (Fig.  5.1). 


A  = 


Supersonic 

Inflow 


Moving  Shock 
-  10  — 


Unsteady 

Subsonic 

Outflow 


Figure  5.1  Unsteady  Diverging  Nozzle  Conhguration 


A  common  steady-state  quasi-l-D  nozzle  test  case  (36)  was  modihed  to  gen¬ 
erate  an  unsteady  flowheld.  Assumed  conditions  at  the  inlet  were:  Mach  number  = 
1.5,  density  =  0.002241  slugs/ft^,  pressure  =  20001b//ft^,  and  temperature  =  520R. 
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Unsteadiness  was  introduced  by  adding  a  sinusoidal  variation  to  the  exit  density. 


p  =  0.003954  +  0.0001  sin(0.l7r  t)  slugs/ft^  (5.7) 

Density  was  varied  as  a  matter  of  convenience-the  desired  unsteady  behavior  was 
produced,  and  a  specihed-density  nozzle  boundary  condition  was  available.  For 
computation,  all  values  were  nondimensionalized  by  the  inlet  velocity  and  density. 

Without  the  unsteady  term  in  the  exit  density,  these  boundary  conditions  result 
in  a  flow  held  with  a  normal  shock  near  the  midpoint  of  the  nozzle  and  subsonic  exit 
how  (36).  With  the  addition  of  the  specihed  unsteadiness  at  the  exit,  the  how  behind 
the  shock  becomes  smoothly  unsteady,  but  remains  subsonic.  The  location  of  the 
shock  oscillates  about  its  steady-state  location,  while  the  supersonic  how  ahead  of 
the  upstream  limit  of  shock  motion  is  not  ahected  and  remains  steady. 

The  unsteady  amplitude  was  selected  by  hrst  hnding  the  range  of  constant 
densities  for  which  a  steady-state  solution  could  be  obtained  using  the  initial  con¬ 
ditions  described  below.  The  amplitude  was  set  to  the  largest  value  that  produced 
upper  and  lower  extremes  within  that  range.  The  frequency  a;  =  O.Itt  was  selected  to 
maximize  the  extent  of  shock  motion  and  still  produce  an  interesting  subsonic  how 
behind  the  moving  shock.  Higher  disturbance  frequencies  produced  a  more  complex 
subsonic  how,  but  resulted  in  a  restricted  range  of  shock  motion.  Lower  disturbance 
frequencies  resulted  in  essentially  the  same  range  of  shock  motion,  but  produced  less 
complex  how  behind  the  shock. 

Flow  solutions  were  obtained  for  three  computational  grids  with  densities  of 
256,  512,  and  1024  cells.  A  uniform  cell  size  was  maintained  in  each  grid.  Solutions 
were  initialized  to  inlet  conditions  over  the  hrst  28%  of  the  grids.  Over  the  remainder 
of  the  grids,  density  and  total  energy  were  initialized  to  inlet  conditions,  while  veloc¬ 
ity  was  initialized  to  34%  of  inlet  velocity.  Adapted  harmonic  balance  solutions  to 
the  unsteady  problem  were  computed  for  six  augmentation  thresholds  from  5.0e-2  to 
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5.0e-7  in  order-of-magnitude  increments.  In  addition,  one  steady-state  (0-frequency) 
solution  was  computed  on  each  grid  to  validate  the  basic  solver. 


Augmentation  increments  available  to  the  adaptation  algorithm  were  specihed 
in  terms  of  candidate  numbers  of  frequencies  that  could  be  included  in  the  solution 
in  a  given  cell.  The  available  candidates  are  listed  in  Table  5.1.  Since  the  flow  in 
much  of  the  grid  was  steady-state,  each  of  the  solutions  was  begun  with  0  frequen¬ 
cies.  The  unsteady  flow  downstream  of  the  shock  was  smoothly  varying,  and  thus 
had  low  frequency  content,  so  all  of  the  lower-frequency  candidates  were  included. 
Finally,  higher  frequency  candidates  in  evenly  spaced  increments  were  included  to 
allow  refinement  of  the  moving  shock  region. 


Table  5.1  Candidate  Numbers  of  Frequencies  Available  to  the  Adaptive  Harmonic 
Balance  Solver  for  the  Unsteady  Nozzle  Test  Case 


0 

1 

2 

3 

4 

5 

6 

7 

10 

13 

16 

19 

22 

24 

27 

32 

37 

42 

A  subsonic  adaptation  scheduling  strategy  was  adopted  (See  Section  4.4).  The 
primary  adaptation  triggers  for  the  initial  and  following  adaptations  were  set  to  2.6 
and  0.1,  respectively,  while  the  secondary,  iteration  based  triggers  were  set  to  7000 
and  1200.  In  addition,  fringe  augmentation  widths  of  4,  8,  and  25  were  required  for 
the  256,  512,  and  1024  cell  grids,  respectively. 

All  solutions  were  initiated  with  a  CFL  of  1.7.  One  test  case  (1024-cell  grid 
with  augmentation  threshold  of  5.0e-4)  required  a  scaled  CFL  to  maintain  stability. 
For  this  case,  the  CFL  was  scaled  down  to  1.25  at  the  highest  frequency  content. 
All  harmonic  balance  solutions  were  converged  to  a  residual  below  l.Oe-6,  a  drop  of 
approximately  4.6  to  4.9  orders  of  magnitude. 

For  comparison  purposes,  conventional  time-accurate  unsteady  solutions  were 
also  computed  on  each  grid.  Each  solution  was  begun  from  a  steady-state  calculation 
based  on  the  mean  outflow  properties,  and  run  time-accurately  for  20  outflow  cycles. 
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Figure  5.2  Time-accurate  Development  of  the  Fully-developed  Unsteady  Nozzle 
Flow  Field 

The  mass  in  the  nozzle  was  tracked  as  the  time-accurate  solution  developed  to  verify 
that  a  fully-developed  solution  was  reached.  Based  on  this  measure,  fully  developed 
flow  was  reached  at  about  the  12th  cycle  (Fig.  5.2). 

5-4  Results 

5.4-1  Comparison  with  Exact  Solutions.  To  validate  the  basic  steady- 
state  quasi- 1-D  harmonic  balance  solver,  0- frequency  (steady-state)  solutions  were 
calculated  for  each  grid  density  and  compared  with  an  exact  solution.  The  exact 
solution  for  a  diverging  nozzle  with  constant  outflow  properties  is  easily  obtained 
from  the  isentropic  flow  equations  and  Rankine-Hugoniot  normal  shock  relations 
(37,  36). 

All  solutions  showed  good  agreement  with  the  theoretical  solution,  as  illus¬ 
trated  in  Fig.  5.3.  The  root-mean-square  (RMS)  percent  error  was  4.9%,  3.5%,  and 
2.5%  for  the  256  cell,  512  cell  and  1024  cell  solutions,  respectively.  Most  of  the 
observed  error  was  attributed  to  smearing  of  the  normal  shock  across  several  cells. 


5-6 


Figure  5.3  Comparison  of  0-Frequency  (Steady  State)  and  Theoretical  Pressures 
for  1024  Dell  Grid  Diverging  Nozzle 

Excluding  the  error  in  the  shock  region  reduced  the  RMS  percent  error  for  the  coarse, 
medium,  and  hne  grids  to  0.31%,  0.16%,  and  0.08%. 

The  validity  of  the  unsteady  harmonic  balance  solutions  was  conhrmed  by 
examining  the  time-average  mass  flux  throughout  the  grid.  With  constant  supersonic 
flow  at  the  inlet  boundary,  conservation  of  mass  requires  that  the  time-averaged  mass 
flux  at  any  given  point  in  the  nozzle  be  inversely  proportional  to  the  ratio  of  nozzle 
cross-sectional  area  at  that  point  to  the  area  at  the  inlet.  In  Fig.  5.4,  the  computed 
time-average  mass  flux  (the  zero-frequency  Fourier  coefficient  of  pu)  for  the  256  and 
1024  cell  grids  is  plotted  along  with  the  theoretical  value.  The  harmonic  balance 
solutions  show  good  agreement  with  the  theoretical  values  everywhere  except  in  the 
region  of  shock  motion,  where  the  effect  of  shock  smearing  introduces  errors.  The 
error  distribution  of  the  time-average  momentum  for  the  256,  512,  and  1024  cell 
harmonic  balance  calculations  is  shown  in  Fig.  5.5.  In  all  three  cases,  the  error 
peaks  at  the  edges  of  the  region  of  shock  motion,  with  the  maximum  error  ranging 
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from  3.5%  for  the  256  cell  solution  to  1.5%  for  the  1024  cell  solution.  Errors  of  this 
magnitude  are  acceptable  for  an  unsteady  analysis  of  this  type. 

5.4-2  Unsteady  Nozzle  Solutions.  The  unsteady-outflow  nozzle  conhgura- 
tion  proved  to  be  a  challenging  test  of  the  adaptive  harmonic  balance  solver.  The 
transition  from  steady-state  flow  to  high-frequency  content  flow  presented  difficul¬ 
ties  during  the  early  transient  stages  of  solution  development.  As  the  solution  was 
rehned,  the  upstream  limit  of  shock  motion  shifted  upstream  slightly.  Higher  frequen¬ 
cies  were  needed  in  the  transition  cells  to  allow  the  disturbance  to  shift  upstream,  but 
none  were  available  to  the  solntion.  Because  the  upstream  flow  was  supersonic  and 
steady-state,  the  transition  cells  had  zero  energy  in  the  hrst  Fourier  frequency,  and 
so  no  freqnencies  were  added  by  previous  threshold-based  augmentation.  In  many 
cases,  the  solntion  in  the  frequency-dehcient  cells  broke  down  before  an  adaptation 
was  triggered  and  more  freqnencies  were  added. 

This  sitnation  was  snccessfnlly  resolved  by  including  fringe  angmentation  (Sec¬ 
tion  2.5).  This  effectively  pre-augmented  a  nnmber  of  cells  ahead  of  the  developing 
nnsteadiness  and  gave  the  developing  unsteadiness  room  to  propagate  upstream. 
(Two  alternative  solutions  to  the  problem — starting  the  solution  with  a  larger  num¬ 
ber  of  frequencies,  and/or  decreasing  the  amonnt  of  flow  development  between  adap¬ 
tations  by  lowering  the  repeat  adaptation  triggers — were  rejected  as  having  too  neg¬ 
ative  an  impact  on  run-time  performance.)  Fringe  augmentation  provided  a  reliable 
solntion  to  the  problem  with  minimnm  increase  in  rnn  time. 

A  typical  growth  history  of  an  adapted  frequency  distribution  is  shown  in 
Fig.  5.6.  The  hgure  shows  the  complete  adaptation  history  for  a  solution  on  the 
hnest  grid  with  an  augmentation  threshold  of  5.0e-7.  Each  frame  shows  the  frequency 
distribution  after  frequency  augmentation  at  the  indicated  iteration,  from  the  hrst 
adaptation  at  iteration  4,725,  to  the  hnal  adaptation  at  iteration  16,990.  As  can  be 
seen  in  Fig.  5.6,  the  frequency  distribution  had  stabilized  ontside  the  shock  region  by 
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Figure  5.6  Evolution  of  the  Adapted  Frequency  Distribution  for  Each  Adaptation 
of  the  Solution,  1024  Cells,  Augmentation  Threshold  5.0e-7 
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the  fourth  adaptation.  Subsequent  adaptations  added  frequencies  only  to  the  region 
of  shock  motion. 

In  Fig.  5.7,  pressure  distributions  for  this  same  adapted  solution  are  shown  at 
10  snapshots  in  time  spanning  one  period  of  the  flow  oscillation.  Also  shown  are 
snapshots  from  a  time-accurate  calculation  on  the  same  computational  grid,  along 
with  the  hnal  frequency  envelope.  This  hgure  illustrates  how  the  adapted  frequency 
distribution  relates  to  the  computed  flow  solution.  Figure  5.7  also  illustrates  how  the 
adapted  harmonic  balance  solution  sharply  captured  the  oscillating  normal  shock 
and  agreed  with  the  time-accurate  solution  at  all  points  in  time.  The  spatially 
varying  frequency  content  produced  no  visible  artifacts  in  the  reconstructed  harmonic 
balance  solution. 

A  quantitative  assessment  of  the  adapted  harmonic  balance  solution  was  ob¬ 
tained  by  calculating  the  RMS  of  the  percentage  differences  between  the  adapted 
solution  and  the  time-accurate  solution  for  each  of  the  10  samples  in  Fig.  5.7,  and 
averaging  the  results.  The  average  RMS  difference  was  0.14%.  (The  maximum  RMS 
difference  over  the  10  samples  was  0.27%.)  A  similar  calculation  was  performed  for 
solutions  obtained  on  all  three  grids  with  all  six  augmentation  thresholds.  The  re¬ 
sults  are  plotted  in  Fig.  5.8.  While  there  was  some  difference  in  the  average  RMS 
values  for  large  augmentation  thresholds,  the  solutions  on  all  three  grids  converged 
to  provide  nearly  identical  agreement  at  a  threshold  of  5.0e-7. 

Nearly  all  of  the  difference  between  the  time-accurate  solutions  and  adapted 
harmonic  balance  solutions  with  high  augmentation  thresholds  was  in  the  region  of 
the  moving  shock.  Figure  5.9  illustrates  how  this  region  developed  with  decreases  in 
augmentation  threshold.  This  hgure  shows  solutions  in  the  region  of  shock  motion 
for  the  256  cell  and  1024  cell  grids  at  a  time  corresponding  to  1/4  of  the  oscillation 
period. 

For  a  given  augmentation  threshold,  the  solutions  on  both  grids  had  similar 
properties.  At  a  threshold  of  5.0e-2,  the  shock  was  poorly  dehned,  and  there  was 
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Figure  5.8  Root-Mean-Square  of  Difference  Between  Adaptive  Harmonic  Balance 
and  Time-accurate  Solutions  with  Changing  Adaptation  Threshold 


Figure  5.9  Time  Snapshot  of  Adaptive  Harmonic  Balance  Solution  Illustrating  Re- 
hnement  of  Oscillating  Normal  Shock  with  Decreased  Augmentation 
Threshold 
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significant  overshoot  behind  the  shock.  At  a  threshold  of  5.0e-4,  the  shock  was  well 
dehned,  but  there  was  ringing  before  and  after  the  shock.  With  a  threshold  of  5.0e- 
6,  the  shock  was  crisply  dehned  on  both  grids,  and  most  of  the  ringing  around  the 
shock  was  removed.  While  the  defects  at  each  threshold  level  were  more  pronounced 
on  the  hner  grid  due  to  higher  grid  resolution  and  frequency  content,  the  basic  form 
of  the  defects  was  the  same  on  both  grids. 

The  frequency  distributions  for  the  solutions  shown  in  Fig.  5.9  are  shown  in 
Fig.  5.10.  These  plots  illustrate  the  relative  number  of  frequencies  required  to  achieve 
each  level  of  solution  hdelity.  They  also  illustrate  how  the  adapted  frequency  distri¬ 
bution  automatically  adjusted  for  grid  density.  Though  a  given  threshold  achieved 
a  qualitatively  similar  solution  on  each  grid,  the  number  of  frequencies  included  in 
the  shock  region  on  the  hue  grid  was  signihcantly  higher  than  the  number  included 
on  the  coarse  grid. 

Figure  5.10  also  includes  the  hnal  frequency  distribution  for  solutions  with 
a  threshold  of  5.0e-7.  As  can  be  seen,  a  large  increase  in  frequency  content  was 
needed  to  achieve  a  small  improvement  in  the  solution  as  indicated  in  Fig.  5.8. 
The  variation  in  maximum  frequency  content  with  grid  density  and  augmentation 
threshold  is  shown  for  all  grids  and  thresholds  in  Fig.  5.11. 

While  it  was  the  maximum  frequency  content  that  determined  the  fidelity  of 
each  solution,  it  was  the  average  frequency  content  that  determined  the  computa¬ 
tional  cost  of  each  solution.  The  average  frequency  content  for  each  adapted  solution 
is  shown  in  Fig.  5.12.  Since  the  discontinuous  flow  was  limited  to  a  small  portion 
of  the  computational  domain,  the  overall  average  frequency  content  remained  small 
for  all  augmentation  thresholds.  On  the  hnest  grid,  average  frequency  content  at 
the  smallest  threshold  grew  by  approximately  4.5  frequencies,  while  the  maximum 
frequency  grew  by  35  frequencies. 

Low  average  frequency  content  translated  directly  into  reduced  compute  times, 
as  shown  in  Fig.  5.13.  As  can  be  seen  in  this  hgure,  the  compute  time  grew  more 
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Figure  5.10  Adapted  Frequency  Distributions  for  Different  Augmentation  Thresh¬ 
olds  on  Coarse  Grid  (top)  and  Fine  Grid  (bottom) 
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Figure  5.11  Change  in  Maximum  Frequency  Content  of  Adapted  Harmonic  Bal¬ 
ance  Solutions  with  Changes  in  Augmentation  Threshold 

slowly  than  the  average  frequency  content.  This  was  due  to  the  inherent  efficiency 
of  resolving  the  flow  features  that  are  slowest  to  converge  during  the  early  stages  of 
the  adaptive  solution  when  the  frequency  content  was  low.  On  the  coarsest  grid,  the 
compute  time  for  the  lowest  augmentation  threshold  is  less  than  twice  that  of  the 
highest  threshold.  On  the  hnest  grid,  the  compute  time  grew  only  2.3  times  from 
the  highest  threshold  to  the  lowest. 

To  determine  the  efficiency  of  the  adaptive  harmonic  balance  approach  relative 
to  non-adapted  harmonic  balance,  each  of  the  harmonic  balance  solutions  was  re¬ 
computed  with  a  hxed  number  of  frequencies.  The  maximum  frequency  content  for 
each  of  the  adapted  solutions  was  identified  and  used  to  generate  the  fixed-frequency 
solutions,  ensuring  that  both  solutions  maintained  the  same  fidelity.  The  efficiency  of 
the  adaptive  method  was  examined  by  comparing  problem  sizes  and  compute  times. 

The  problem  size  for  the  fixed  frequency  calculation  was  determined  by  the 
number  of  frequencies  included  in  the  solution.  The  equivalent  problem  size  for  the 
adapted  calculation  is  given  by  the  average  frequency  content.  The  difference  be- 
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tween  these  two  quantities  is  shown  in  Figure  5.14  in  terms  of  percent  reduction. 
For  the  largest  augmentation  threshold,  the  reduction  in  problem  size  was  approxi¬ 
mately  45%  for  all  grid  densities.  This  represents  a  signihcant  reduction,  considering 
that  the  average  number  of  frequencies  in  each  solution  was  approximately  one,  the 
minimum  possible  number  of  frequencies  for  a  frequency  augmentation  approach. 
Any  improvement  would  require  a  frequency  decimation  approach  where  frequencies 
are  identihed  as  unnecessary  and  removed  from  the  solution.  Such  an  approach  was 
not  implemented  in  the  current  research.  The  problem  size  reduction  was  more  sub¬ 
stantial  for  the  small  thresholds,  approaching  90%  for  the  smallest  threshold  on  the 
hnest  grid. 


Figure  5.14  Reduction  in  Average  Frequency  Content  for  Adapted  vs  Non-adapted 
Harmonic  Balance  Solutions 

The  difference  between  compute  times  for  the  adapted  and  hxed-frequency 
solutions  is  shown  in  Fig.  5.15,  again  in  terms  of  percent  reduction.  It  is  clear  from 
these  results  that  a  given  reduction  in  problem  size  did  not  necessarily  result  in  an 
equivalent  reduction  in  compute  time,  especially  for  higher  augmentation  thresholds 
and  lower  grid  densities.  Dehning  an  efficiency  of  the  adaptive  solver  as  the  ratio 
of  the  compute  time  reduction  to  the  problem  size  reduction,  the  average  efficiency 
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Figure  5.15  Reduction  in  Compute  Time  for  Adapted  vs  Non-adapted  Harmonic 
Balance  Solutions 

for  a  threshold  of  5.0e-2  was  only  47%.  For  lower  thresholds  the  efficiency  was  much 
higher,  averaging  91%  for  a  threshold  of  5.0e-6,  and  95%  for  a  threshold  of  5.0e-7. 
The  adaptive  solution  on  the  hnest  grid  with  an  augmentation  threshold  of  5.0e-7 
achieved  a  99%  adaptive  efficiency  with  an  86%  reduction  in  compute  time.  The 
adaptive  efficiencies  for  all  solutions  are  shown  in  Fig.  5.16. 

5. 5  Summary 

The  adaptive  split-domain  harmonic  balance  method  developed  and  analyzed 
in  previous  chapters  was  successfully  applied  to  solve  the  quasi-l-D  inviscid  flow  in  a 
supersonic-subsonic  diverging  nozzle  with  unsteady  outflow  properties.  This  problem 
is  representative  of  the  types  of  problems  the  adaptive  harmonic  balance  approach 
was  designed  for:  a  stationary  unsteady  flow  held  containing  mostly  smooth,  low- 
frequency  content  how,  but  with  isolated  regions  of  highly  nonlinear  high-frequency 
content  how. 
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Figure  5.16  Efficiency  of  Adaptive  Harmonic  Balance  Solutions  (Run  Time  Re¬ 
duction/Average  Frequency  Reduction) 

Adapted  harmonic  balance  solutions  were  generated  for  a  variety  of  grid  den¬ 
sities  and  augmentation  thresholds.  The  energy-based  frequency  augmentation  ap¬ 
proach  proved  effective  in  matching  frequency  content  to  underlying  flow.  Regions 
of  the  computational  domain  requiring  higher  frequency  content  were  identihed,  and 
frequencies  were  added  until  augmentation  threshold  levels  were  achieved.  Given 
thresholds  produced  qualitatively  uniform  solutions  across  tested  grid  densities.  At 
the  lowest  augmentation  thresholds,  the  adapted  harmonic  balance  solutions  showed 
good  agreement  with  both  theoretical  and  time-accurate  numerical  solutions. 

The  combination  of  an  energy-based  frequency  augmentation  approach,  fringe 
augmentation,  and  the  split-domain  solver  proved  to  be  robust  and  stable.  The 
adaptation  algorithm  and  solution  scheme  successfully  handled  a  rapid  transition 
from  steady-state  supersonic  flow  to  the  highly  nonlinear,  unsteady,  mixed  super¬ 
sonic/subsonic  flow  in  the  region  of  the  moving  shock.  In  all  but  one  case,  the 
maximum  stable  CFL  of  the  solver  was  able  to  be  used  for  the  entire  solution. 
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The  adaptive  harmonic  balance  approach  effectively  reduced  the  time  needed 
to  obtain  a  high-hdelity  solution  with  the  harmonic  balance  method.  The  compute 
time  for  solutions  with  thresholds  of  5.0e-6  and  5.0e-7  averaged  greater  than  80% 
reductions  relative  to  equivalent  non-adapted  solutions  on  all  grids.  With  these 
reductions,  the  difference  in  compute  time  between  low-hdelity  and  high-hdelity 
solutions  was  also  reduced  so  that  the  highest  hdelity  solutions  took  only  2-2.5 
times  longer  than  the  lowest-hdelity  solutions. 


5-21 


VI.  Summary  and  Conclusions 


6. 1  Summary  of  Research 

A  new  adaptive  split-domain  harmonic  balance  CFD  method  was  developed 
and  applied  to  a  variety  of  one-dimensional  problems.  The  new  method  employed 
a  unique  multi-domain  split-operator  (split-domain)  solution  approach  to  efficiently 
solve  harmonic  balance  equations  with  large  numbers  of  Fourier  frequencies.  The 
split-domain  approach  successfully  removed  a  numerical  stability  restriction  present 
in  previous  harmonic  balance  CFD  implementations.  It  also  reduced  the  computa¬ 
tional  work  required  to  calculate  a  solution  by  reducing  the  number  of  FFTs  required 
to  just  one  per  point,  per  iteration.  Improved  stability  and  reduced  computational 
work  both  resulted  in  reduced  run  time. 

To  further  improve  performance,  the  split-domain  solution  method  was  com¬ 
bined  with  a  frequency-adaptation  algorithm  that  minimized  the  size  of  the  harmonic 
balance  problem  by  varying  the  number  of  frequencies  included  in  the  harmonic  bal¬ 
ance  solution.  A  frequency-augmentation  adaptation  approach  was  employed.  With 
this  approach,  harmonic  balance  solutions  were  begun  with  a  minimum  number  of 
frequencies,  and  periodically  examined  to  determine  where  additional  frequencies 
were  needed  to  capture  the  local  flow.  The  decision  to  add  frequencies  was  based  on 
an  examination  of  the  fraction  of  spectral  energy  contained  in  the  largest  included 
Fourier  frequency.  Frequencies  were  added  to  cells  that  contained  more  energy  in 
the  highest  frequency  than  a  user-specihed  maximum  threshold.  By  changing  the 
threshold,  the  number  of  frequencies  in  the  solution,  and  thus  the  solution  accuracy 
and  run  time,  were  controlled. 

The  non-adaptive  split-domain  harmonic  balance  CFD  method  was  applied 
to  the  1-D  inviscid  Burgers’  equation.  Large  amplitude,  time-periodic  solutions 
were  computed  and  compared  to  solutions  computed  with  prior  harmonic  balance 
approaches.  The  split-domain  method  produced  solutions  comparable  to  those  of 
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the  prior  methods,  while  successfully  eliminating  a  stability  restriction  experienced 
by  those  methods  when  a  large  number  of  Fourier  frequencies  are  included  in  the 
solution.  It  was  found  that  the  difference  between  harmonic  balance  solutions  and 
conventional  time-accnrate  solntions  was  asymptotic  with  respect  to  the  number  of 
Fourier  frequencies  included  in  the  harmonic  balance  solution.  When  the  number  of 
frequencies  was  equal  to  or  greater  than  some  asymptotic  frequency,  the  harmonic 
balance  solntions  were  comparable  to  the  time-accnrate  solntions.  Several  factors 
were  found  to  inflnence  the  asymptotic  freqnency,  including  disturbance  frequency, 
the  strength  of  the  moving  discontinnity,  and  the  computational  grid  density. 

The  fnll  adaptive  split-domain  harmonic  balance  CFD  method  was  applied  to 
the  1-D  Euler  equation.  The  method  was  snccessfnlly  employed  to  solve  a  variety  of 
supersonic  and  snbsonic  one-dimensional  flow  fields  containing  strong  moving  shocks. 
The  energy-based  angmentation  approach  reliably  identihed  cells  where  additional 
freqnency  content  was  needed  and  could  be  supported  by  the  computational  grid. 
The  resulting  adapted  harmonic  balance  solutions  were  eqnivalent  to  non-adapted 
harmonic  balance  solutions,  and  compared  well  with  conventional  time-accnrate  so¬ 
lutions. 

Adaptation  schednling  was  found  to  have  a  signihcant  impact  on  the  rnn-time 
performance  of  the  adaptive  solver.  In  some  cases,  applying  the  wrong  scheduling 
approach  more  than  doubled  solution  compute  time.  For  supersonic  flows,  rapid 
adaptation  with  minimal  flow  development  time  produced  the  best  performance.  For 
subsonic  flows,  a  scheduling  approach  that  allowed  more  time  for  flow  development 
between  adaptations  was  best. 

The  compatibility  of  the  adaptive  split-domain  harmonic  balance  approach  and 
FAS  mnltigrid  acceleration  was  demonstrated.  Solutions  were  computed  with  and 
withont  frequency  angmentation,  and  with  and  withont  mnltigrid  acceleration.  The 
relative  performance  benehts  of  each  approach  remained  consistent. 
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Finally,  the  adaptive  split-domain  approach  was  demonstrated  by  solving  for 
quasi-l-D  flow  in  a  supersonic-subsonic  diverging  nozzle  with  unsteady  outflow  prop¬ 
erties.  This  problem  contained  a  wide  range  of  flow  regimes,  including  supersonic 
steady-state  flow,  continuous,  unsteady  subsonic  flow,  and  an  oscillating  normal 
shock.  Adapted  harmonic  balance  solutions  were  generated  for  a  variety  of  grid 
densities  and  augmentation  thresholds.  The  energy-based  frequency  augmentation 
approach  proved  effective  in  matching  frequency  content  to  the  underlying  flow. 
Regions  of  the  computational  domain  requiring  higher  frequency  content  were  iden- 
tihed,  and  frequencies  were  added  until  augmentation  threshold  levels  were  achieved. 
However,  it  was  necessary  to  supplement  threshold-based  augmentation  with  fringe 
augmentation  in  order  to  successfully  handle  shifting  frequency  transitions  during 
the  early  stages  of  solution  development.  Given  thresholds  produced  qualitatively 
uniform  solutions  across  tested  grid  densities.  At  the  lowest  augmentation  thresh¬ 
olds,  the  adapted  harmonic  balance  solutions  showed  good  agreement  with  both 
theoretical  and  time-accurate  numerical  solutions. 

The  adaptive  harmonic  balance  approach  effectively  reduced  the  time  needed 
to  obtain  an  accurate  solution  with  the  harmonic  balance  method.  The  compute 
time  for  solutions  with  thresholds  of  5.0e-6  and  5.0e-7  averaged  greater  than  80% 
reductions  relative  to  equivalent  non-adapted  split-domain  solutions  on  all  grids. 
With  these  reductions,  the  difference  in  compute  time  between  low-hdelity  and  high- 
hdelity  solutions  was  also  reduced,  to  the  point  where  the  highest  hdelity  solutions 
took  only  2-2.5  times  longer  than  the  lowest-hdelity  solutions. 

6.2  Conclusions 

As  a  result  of  the  research  conducted  and  documented  in  this  dissertation,  the 
following  major  conclusions  were  reached: 

1.  The  harmonic  balance  method  produces  accurate  solutions  to  strongly  nonlinear 
time-periodic  flow  problems,  provided  that  a  sufficient  number  of  Fourier  fre- 
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quencies  are  included.  There  is  a  practical  asymptotic  limit  that  is  dependent 
on  the  flow  fleld  and  on  the  density  of  the  computational  grid.  Thus,  adding 
frequencies  does  not  automatically  improve  a  harmonic  balance  solution. 

2.  A  split- operator,  multi-domain  (split- domain)  solution  approach  can  he  em¬ 
ployed  to  efficiently  solve  a  harmonic-balance  CFD  problem  with  a  sufficient 
number  of  frequencies  to  accurately  model  a  strongly  nonlinear  periodic  flow. 
The  split-domain  approach  effectively  reduces  or  removes  stability  restrictions 
associated  with  high-frequency  harmonic  balance  solutions.  It  also  reduces  the 
required  number  of  Fourier  transform  pairs  to  one,  regardless  of  the  time  in¬ 
tegration  scheme  implemented.  Finally,  the  approach  is  fully  compatible  with 
FAS  multigrid  convergence  acceleration. 

3.  The  number  of  frequencies  included  in  a  split-domain  harmonic  balance  solu¬ 
tion  can  be  varied  from  cell  to  cell  without  negatively  impacting  solution  qual¬ 
ity.  With  proper  implementation,  transitions  between  regions  with  different 
frequency  content  are  transparent  and  do  not  affect  a  reconstructed  solution. 
Frequency  transitions  cannot  be  treated  as  time-lagged  boundaries;  both  sides 
of  the  transition  must  be  synchronously  integrated  in  time.  Resampling  via 
truncation  or  zero-padding  of  Fourier  coefficients  was  an  effective  way  to  com¬ 
pute  differences  across  frequency  transitions. 

4.  For  flows  with  time-responses  that  are  mostly  continuous  with  a  finite  number 
of  discontinuities  and  no  impulses,  the  fraction  of  spectral  energy  contained  in 
the  highest  computed  Fourier  frequency  of  a  harmonic  balance  solution  (En) 
is  an  accurate  indicator  of  how  completely  the  solution  has  modeled  local  flow 
behavior.  Solutions  containing  an  insufficient  number  of  Fourier  frequencies 
will  have  a  higher  Em  than  those  that  contain  a  sufficient  number.  When  a 
computational  grid  is  too  coarse  to  support  high  frequencies,  the  Em  is  reduced. 
A  given  minimum  solution  fidelity  can  be  maintained  by  selecting  an  upper 
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bound  on  Ej^j  and  including  sufficient  frequencies  to  ensure  is  below  this 
bound. 

5.  A  spatially- adaptive  harmonic  balance  method  can  he  implemented  to  accurately 
and  efficiently  compute  a  stationary  time-periodic  flow  field  containing  regions 
of  smooth  flow  and  regions  with  strong  moving  discontinuities  by  automatically 
adjusting  the  number  of  freguencies  in  the  solution,  on  a  point-by-point  (cell- 
by-cell)  basis,  to  match  local  flow  conditions.  A  restatement  of  the  research 
thesis.  An  adaptive  split-domain  harmonic  balance  CFD  solver  employing  an 
energy-based  frequency  augmentation  adaptation  approach  was  demonstrated 
and  shown  to  signihcantly  reduce  the  time  required  to  compute  an  accurate 
harmonic  balance  solution. 
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VII.  Recommendations  for  Future  Research 

The  next  logical  step  in  the  development  of  the  adaptive  split-domain  harmonic  bal¬ 
ance  approach  is  the  extension  of  the  method  to  higher  dimensions  and  application 
to  turbomachinery  problems.  Extension  of  the  basic  approach  should  be  straightfor¬ 
ward,  requiring  only  minor  modihcations  to  the  current  implementation.  Part  of  this 
effort  should  be  the  development  of  a  boundary  condition  to  model  a  rotor-stator 
interface.  Additional  areas  that  merit  further  research  are  discussed  below. 

Combined  Grid/ Frequency  Adaptation 

Grid  adaptation  provides  a  means  of  efficiently  resolving  flow  features  when 
their  precise  location  is  not  known  a  priori.  Applying  grid  adaptation  to  conventional 
unsteady  CFD  problems  is  difficult,  however,  because  flow  features  are  in  constant 
motion.  Since  the  harmonic  balance  method  solves  a  steady-state  problem,  imple¬ 
mentation  of  grid  adaptation  should  be  greatly  simplihed.  One  of  the  primary  goals 
of  research  into  combined  grid  and  frequency  adaptation  would  be  to  determine  if 
frequency  content  information  could  be  used  along  with  flow  feature  information  to 
drive  adjustments  to  the  computational  grid. 

Improved  Multigrid  for  Harmonic  Balance  CFD 

One  of  the  primary  incentives  for  developing  harmonic  balance  CFD  approaches 
is  that  a  steady-state  problem  is  solved,  which  allows  convergence  acceleration  tech¬ 
niques  such  as  multigrid  to  be  employed.  The  current  research  demonstrated  the 
compatibility  of  FAS  multigrid  and  adaptive  split-domain  harmonic  balance,  but 
identifying  the  best  multigrid  approach  was  beyond  scope.  Identifying  an  optimal 
multigrid  approach  for  harmonic  balance  problems  would  be  of  great  value. 

One  property  of  a  harmonic  balance  solution  that  limits  the  effectiveness  of 
multigrid  is  that  the  real  and  imaginary  parts  of  the  high  frequency  Fourier  coeffi- 
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cients  tend  to  be  highly  oscillatory  in  the  spatial  dimensions.  Conventional  multigrid 
techniques  are  most  effective  when  the  spatial  behavior  of  the  solution  is  smooth, 
so  that  little  information  is  lost  when  the  solution  is  transfered  to  a  coarser  grid. 
When  multigrid  is  applied  to  oscillatory  problems,  information  is  lost  during  grid 
transfer  operations,  and  convergence  acceleration  is  reduced.  Some  attempts  have 
been  made  to  develop  multigrid  techniques  for  oscillatory  elliptic  PDEs  (e.g.  (38)). 
One  avenue  of  research  would  be  to  adapt  these  techniques  to  the  harmonic  balance 
equations. 

Another  approach  suggested  by  the  current  research  would  be  to  transfer  so¬ 
lutions  between  grids  using  the  magnitude-phase  form  of  the  Fourier  coefficients, 
rather  than  the  real  and  imaginary  form.  It  was  observed  that  the  magnitude  of 
the  complex  Fourier  coefficients  was  considerably  less  oscillatory  than  the  real  and 
imaginary  components.  The  phase,  when  “unwrapped”  so  that  it  was  continuous  and 
had  values  greater  than  27r,  was  nearly  linear  with  respect  to  the  spatial  dimension. 
The  smoother  properties  of  a  magnitude/phase  should  allow  more  accurate  transfer 
between  grids.  The  major  challenge  facing  such  an  approach  is  developing  a  reliable 
means  of  unwrapping  phase. 

Magnitude/ Phase  Harmonic  Balance  Solver 

In  order  to  accurately  capture  the  oscillatory  behavior  of  the  real  and  imagi¬ 
nary  coefficients  of  higher-frequency  Fourier  coefficients,  a  hue  computational  grid  is 
necessary.  As  noted  above,  spatial  variation  of  the  magnitude  and  unwrapped  phase 
of  the  coefficients  is  much  smoother,  and  thus  can  be  captured  with  a  much  coarser 
grid  (39).  Development  of  a  harmonic  balance  method  based  on  a  magnitude/phase 
representation  of  the  complex  coefficients,  could  signihcantly  reduce  the  number  of 
grid  cells  required  to  compute  an  accurate  harmonic  balance  solution,  and  thus  re¬ 
duce  the  memory  and  run-time  requirements  as  well.  A  magnitude/phase  solver 
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would  also  improve  multigrid  efficiency,  increasing  the  accuracy  of  the  corrections 
computed  on  coarse  grids. 

Many  aspects  of  a  magnitude/phase  solver  would  be  similar  to  the  current  split- 
domain  solver.  Boundary  condition  and  flux  calculations  could  still  be  done  in  time 
domain.  However,  face  fluxes  would  likely  be  computed  from  averaged  cell  magnitude 
and  phase,  instead  of  averaged  cell  fluxes  as  in  the  current  implementation.  As  with 
magnitude/phase  multigrid  transfers,  a  major  research  challenge  is  determining  how 
to  correctly  implement  continnons  phase. 


Efficient  Single-domain  Harmonic  Balance 

The  snggested  research  areas  discussed  thus  far  have  all  involved  the  extension 
or  modihcation  of  the  adaptive  split-domain  harmonic  balance  method.  An  inves¬ 
tigation  of  an  efficient  single-domain  harmonic  balance  approach,  while  snggested 
by  the  split-domain  approach,  is  a  research  area  that  could  result  in  a  signihcant 
departnre  from  the  cnrrent  method. 

The  basis  of  the  proposed  research  is  the  physical  interpretation  of  the  split- 
domain  solntion  process  discussed  in  Section  2.4.  It  was  shown  that  the  split-domain 
solntion  approach  is  eqnivalent  to  a  sequence  of  snccessive  integrate-shift  operations 
in  which  Eq.  2.23a,  repeated  here  as  Eq.  7.1,  is  integrated  forward  in  psendo  time, 
and  the  elements  of  the  solution  vector  f  are  shifted  backward  in  physical  time. 


dj  a«i-(0 

dr  dx 


(7.1) 


A  steady-state  solution  is  achieved  when  the  change  in  f  dne  to  the  psendo-time 
integration  step  is  the  negative  of  the  change  dne  to  the  physical  time  shift. 


The  combined  size  of  the  physical-time  shifts  was  shown  to  be  equal  to 

At,  where  At  is  the  physical  time  sample  increment.  If  Ar  =  At  =  (^2n+i)oj  ’ 
then  the  size  of  the  physical  time  shift  is  eqnal  to  the  integration  step  size.  At.  In 
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this  case,  the  distinction  between  pseudo-time  and  physical  time  disappears,  and  the 
change  in  ^  due  to  the  shift  operation  is  simply  the  difference  between  consecutive 
elements  of  For  this  integration  step  size,  there  is  no  need  to  transform  to  the 
frequency  domain  to  compute  the  shift. 

Instead  of  simply  fixing  the  time  integration  step  size  and  replacing  the  fre¬ 
quency  domain  operations,  an  alternative  form  of  the  harmonic  balance  equations 
could  be  employed.  If  pseudo-time  is  replaced  by  physical  time  and  Eq.  7.1  is  inte¬ 
grated  from  an  initial  condition  at  time  t  to  time  t  +  At,  then  the  harmonic  balance 
problem  reduces  to 


where  A  is  a  2iV  -h  1  x  2iV  -h  1 


A  = 


This  approach  would  be  more 
required. 

Both  single-domain  approaches  eliminate  the  need  for  Fourier  transforms  dur¬ 
ing  the  solution  process,  and  thus  might  reduce  the  time  required  to  compute  a 
solution.  The  second  approach  has  the  additional  advantage  that  the  frequency  of 
oscillation,  u,  does  not  appear  in  the  equations,  and  thus  would  not  need  to  be 
known  a  priori]  it  would  enter  the  solution  only  through  the  boundary  conditions. 


A 


dx 


=  0, 


matrix  with  the  form 


(7.2) 


1  0  ...  0  -1 

-110  0 

0  -1  1  ; 

;  0 

0  ...  0  -1  1 


(7.3) 


flexible,  because  no  specific  integration  time  step  is 
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Appendix  A.  Derivation  of  Harmonie  Balanee  Burgers’  Equation  for 

a  Complex  Fourier  Series 

In  this  appendix,  a  direct-substitution  implementation  of  the  one-dimensional  har¬ 
monic  balance  Burgers’  equation  is  derived. 

The  one-dimensional  inviscid  Burgers  equation  is  given  by 


du  1  dC 
Ft  ^  2'dx 


(A.l) 


where  u{x,  t)  is  the  dependent  variable,  t  is  the  temporal  variable,  and  x  is  the  spatial 
variable.  Let  u{x,t)  be  approximated  by  a  truncated  complex  Fourier  series 

N 

u{x,t)  ^  Cn(a:)e*"^*  (A.2) 

n=-N 


where  N  is  the  number  of  positive  frequencies  in  the  truncated  series,  and  on  is  the 
fundamental  frequency  of  the  series.  In  this  series,  the  complex  coefficients,  an, 
are  functions  of  the  spatial  variable  only.  Substituting  Eq.  A.2  into  Eq.  A.l  and 
evaluating  the  time  derivative,  one  obtains 


N 


n=-N 


2  ox 


j  =  0. 


(A.3) 


As  a  result  of  squaring  the  approximating  series  in  the  flux  term,  Eq.  A.3  contains 
high  frequency  terms  with  \m  +  n\  >  N .  As  part  of  the  harmonic  balance  approx¬ 
imation,  these  terms  are  discarded,  leaving  only  terms  with  frequencies  included  in 
the  approximating  Eourier  series.  Gathering  the  remaining  terms  of  like  frequency 
gives 


N 

E 

fc=0 


ikojCk  + 


2  dx 


N 


N 


(m+n),/c 


\m=—N  n=—N 


(A.4) 
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In  these  equations,  S(^m+n),k  is  the  Kronecker  delta  function,  dehned  as  1  if  (m  +  n)  = 
k,  and  0  if  (m  +  n)  7^  k,  thus  the  summation  inside  the  spatial  derivative  contains 
only  products  of  Cm  and  c„  such  that  m  +  n  =  k. 

The  frequencies  in  Eq.  A. 4  are  now  “balanced”.  Since  the  right  hand  side  is 
identically  zero,  this  simply  means  that  each  of  the  terms  in  the  outer  summation  is 
required  to  satisfy  equality  individually.  For  the  resulting  2N  +  1  complex  equations 
to  hold  for  all  time,  the  terms  inside  the  square  brackets  must  be  identically  equal 
to  zero.  The  exponential  terms  can  be  dropped,  leaving  a  system  of  2N  +  1  coupled 
ordinary  differential  equations  for  2N  +  1  complex  coefficients.  Since  the  dependent 
variable,  u{x,t),  is  real  and  periodic,  the  negative  frequency  Fourier  coefficients  are 
just  the  complex  conjugates  of  the  corresponding  positive  frequency  coefficients,  i.e. 
C-n  =  Cn-  Furthermore,  the  steady-state  coefficient  cq  must  be  real,  bringing  the  total 
number  of  unknowns  to  2A^-|-1:  the  real  and  imaginary  parts  of  N  positive-frequency 
coefficient,  plus  the  steady-state  coefficient. 

Keeping  just  the  positive  frequency  equations  and  eliminating  the  exponential 
terms,  Eq.  A. 4  reduces  to 


ikojCk  + 


lA 

2  dx 


N 


N 


^  ^  ^  ^  ^ (m+n) ,k  0 


0  <  A;  <  A^ 


(A.5) 


\m=—N  n=—N 


or  in  vector  form. 


where 


S(U)  + 


dF(U) 

dx 


=  0 


(A.6) 


Co 

0 

2^m=-N  2^n=-N  C(m+n),0 

u  = 

Cl 

S(U)  = 

iuci 

F(U)  =  i 

Z^m=-N  2^n=-N  C(m+n),l 

Cn 

iNujcn 

2^m=-N  l^n=-N  ^(m+n),N 

(A.7) 
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and  is  defined  as  c„  if  n  >  0,  and  c„  if  n  <  0.  An  example  of  the  harmonic  balance 
Burgers’  equation  flux  vector  is  presented  here  for  N  =  3. 


F(U)^=3 


|(Co  +  |cip  +  |c2p  +  |c3p) 

CqCi  +  C1C2  +  C2C3 
C0C2  +  C1C3 
C0C3 


(A.8) 
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Appendix  B.  Impact  of  Operator  Splitting  Error  in  the  Split-Domain 

Harmonic  Balance  Solution 

Introduction 

One  potential  problem  with  a  split-operator  solution  technique  such  as  the 
split-domain  harmonic  balance  formulation,  is  that  for  nonlinear  problems,  the  op¬ 
erator  splitting  is  not  exact  (23).  In  addition,  discretization  of  a  split  equation 
typically  results  in  the  introduction  of  an  error  that  is  dependent  on  the  size  of  the 
numerical  time  step  taken  to  solve  the  problem.  Therefore,  it  is  necessary  to  under¬ 
stand  the  nature  of  the  splitting  error  for  a  specihc  application  before  the  approach 
is  used.  This  appendix  examines  the  error  introduced  by  the  split-domain  harmonic 
balance  approach  both  analytically  and  experimentally. 


Linear  analysis  of  Splitting-induced  Error 
Consider  the  linear  scalar  equation 

+  afx  —  bf  =  0 


(B.l) 


where  f  and  b  are  complex  scalars,  and  fx  is  the  spatial  derivative  of  f.  Splitting 
into  two  homogeneous  equations  gives 


(B.2) 

(B.3) 


For  simplicity,  assume  a  first-order  forward-Euler  discretization  in  time.  The 
effective  hnite  difference  equation  resulting  from  a  symmetric  Strang  splitting  which 
evaluates  the  ODE  twice  is  obtained  by  successively  applying  the  time  discretization. 
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Integrate  the  ODE  1/2  time  step: 


r 


(B.4) 


Integrate  the  PDE  a  full  time  step  using  the  results  of  the  previous  integration: 


r*  =  r-A«(o(r)j 

=  ^  -  AioK  + 

b  /\f^ 

=  i  +  -  aCx) - ^ab^x  (B.5) 


Finally,  integrate  the  ODE  an  additional  1/2  time  step,  again  using  the  results  of 
the  previous  integration  step: 


*n+l  _ 


r-  + 


A  2  fb^  ^  a6 

-  ire 


At_ 

^o-b^x  +  ^^e  + 
At^ 


=  e  +  At(6e  -  a^x)  +  At^  -e  -  -  —ab^i 


ab'^ix 
62 


At^ 
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Subtracting  ^  +  At{b^  —  a^x)  from  both  sides  of  Eq.  B.6  and  dividing  by  At 


gives 


cn+l  _  f  /la  \  \f2 

— +  aix  -bi  =  At  \  —i-  abix  j  -  —ab'^^x 


(B.7) 


At  steady-state,  the  hrst  term  on  the  left  hand  side  is  equal  to  zero.  The  semi- 
discretized  equation  satished  at  steady-state  is  therefore  given  by 


^e*  ^e 


At 


62 


e  -  ab^a 


+  H.O.T. 


(B.8) 
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The  first  term  on  the  right  hand  side  of  Eq.  B.8  represents  the  time-step- 
dependent  splitting  error.  More  insight  into  the  effect  of  this  error  is  gained  by 
discarding  the  higher  order  terms  and  rewriting  the  equation. 

/  hi^  /\t\ 

{bAt  +  l)a^^-  \b+ ^  =  0  (B.9) 

This  shows  that  the  splitting  error  affects  both  the  wave  speed  and  the  strength 
of  the  source  term.  Since  the  current  research  involved  only  explicit  solvers,  bAt  is 
assumed  to  be  much  less  than  a,  and  the  error  in  wave  speed  is  disregarded.  The 
remaining  discussion  will  concentrate  on  the  source  term  error  component. 

The  harmonic  balance  source  term  is  given  by 

0 

..  ..  iuci 

S{0  =  ,  .  (B.IO) 

iNojcn 

The  harmonic  balance  source  term  has  the  same  form  as  the  source  in  the  model 
equation,  with  b  =  inoj.  Assuming  the  results  of  the  simplihed  analysis  can  be 
applied  to  each  term  independently,  the  effective  source  term  for  the  split-domain 
harmonic  balance  equations  becomes 

0 

{-^  +  iu)ci 

The  splitting  error  adds  a  negative  real  component  to  the  harmonic  balance  source 
coefficient.  Several  conclusions  can  be  drawn  from  this  result.  First,  it  is  clear 
that  the  steady-state  (zero-frequency)  equation  is  unaffected  by  splitting  error.  The 


(B.ll) 
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error  should  remain  small  for  low  frequencies,  but  will  increase  rapidly  in  the  high 
frequencies  due  to  the  term.  Since  the  error  term  contains  the  square  of  the 
fundamental  frequency,  the  splitting  error  will  be  somewhat  more  prominent  for 
high-frequency  unsteady-flows. 

Experimental  Analysis  of  Splitting-indueed  Error 

To  assess  the  actual  effect  of  splitting  error  on  a  harmonic  balance  solution, 
supersonic  and  subsonic  solutions  (conhgurations  SS3  and  SUl  from  Chapter  IV) 
were  computed  at  a  variety  of  CFLs.  Each  of  the  solutions  was  computed  on  a  601 
point  grid.  The  supersonic  solutions  were  computed  with  32  frequencies,  while  the 
subsonic  solutions  were  computed  with  45  frequencies. 

Figures  B.1-B.4  show  a  comparison  of  the  real  part  of  the  5th,  15th,  25th,  and 
32nd  Fourier  coefficients  of  momentum,  respectively,  for  the  supersonic  case  at  CFLs 
of  0.5,  1.1,  and  1.7.  As  predicted  by  the  simple  linear  analysis,  there  was  little  dif¬ 
ference  between  the  solutions  for  the  low  frequency  coefficient  (Fig.  B.l).  For  higher 
Fourier  frequencies,  the  effect  of  splitting  error  became  more  pronounced  (Figs.  B.2, 
B.3,  and  B.4).  As  the  splitting  error  increased  with  larger  time  steps  (CFLs),  the 
solution  became  more  dissipative.  The  frequency  of  the  spatial  oscillations  in  each 
coefficient  remained  essentially  unchanged.  Results  for  the  imaginary  part  of  the 
Fourier  coefficients  were  similar. 

Figures  B.5-B.9  show  similar  comparisons  of  the  real  part  of  the  5th,  15th, 
25th,  35nd  and  45th  Fourier  coefficients  of  momentum  for  the  subsonic  case  at  CFLs 
of  0.6,  1.0,  and  1.4.  For  low  to  mid  frequencies  (Figs.  B.5,  B.6,  and  B.7),  the 
splitting  error  exhibits  the  same  increasingly  dissipative  effect  as  in  the  supersonic 
case.  However,  at  higher  frequencies,  the  dissipative  effect  of  the  splitting  error  seems 
to  decrease  again. 
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The  discussion  so  far  has  focused  on  the  impact  of  splitting  error  on  individual 
Fourier  coefficients.  Figures  B.IO  and  B.ll  illustrate  how  the  splitting  error  affects 
the  reconstructed  solutions. 

As  can  be  seen  in  Fig.  B.IO,  splitting  error  had  little  effect  on  the  supersonic 
solution.  Solutions  at  the  lowest  and  highest  CFL  were  qualitatively  similar  except 
for  the  presence  of  increased  oscillations  behind  the  moving  shocks  with  smaller  time 
steps.  The  subsonic  solution  also  experiences  an  increase  in  oscillations  behind  the 
moving  shocks  with  a  smaller  time  step  (Fig.  B.ll).  However,  in  the  trough  re¬ 
gions,  the  low-CFL  solution  shows  a  slight  improvement  over  the  high-CFL  solution, 
partially  overcoming  the  error  introduced  at  the  exit  boundary  (See  Section  4.5). 

For  both  the  supersonic  and  subsonic  cases,  the  oscillations  that  occurred  be¬ 
hind  the  shocks  for  low-CFL  solutions  were  not  readily  removed  by  increasing  the 
amount  of  artificial  dissipation  in  the  solution.  Any  dissipation  strong  enough  to  af¬ 
fect  these  oscillations  also  damped  low-frequency  modes  that  were  not  significantly 
affected  by  changes  in  time  step  size,  resulting  in  an  overall  degradation  of  the  solu¬ 
tion. 

Conclusions 

A  simplihed  analysis  of  splitting-induced  error  in  the  split-domain  harmonic 
balance  method  suggested  that  any  time-step  size  dependency  in  a  converged  solution 
would  be  most  evident  in  the  higher-frequency  terms.  Just  such  an  error,  in  the  form 
of  damping  of  high-frequency  Fourier  coefficients,  was  observed  experimentally  for 
both  supersonic  and  subsonic  flow  helds.  The  observed  error  had  little  effect  on 
reconstructed  flow  properties  such  as  pressure.  If  anything,  the  increased  splitting 
error  resulting  from  a  larger  time  step  had  a  positive  rather  than  negative  effect  on 
the  solution,  eliminating  high-frequency  oscillations  that  were  otherwise  difficult  to 
remove. 
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Figure  B.l 


Effect  of  Splitting  Error  on  5*^  Fourier  Coefficient  with  Changing  CFL 
for  Supersonic  Flow 


vied  leey  ‘luepmeoQ  Lun}U9Luo|/\| 


Figure  B.2 


Effect  of  Splitting  Error  on  15*^  Fourier  Coefficient  with  Changing  CFL 
for  Supersonic  Flow 
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Figure  B.3  Effect  of  Splitting  Error  on  25*^  Fourier  Coefficient  with  Changing  CFL 
for  Snpersonic  Flow 
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Figure  B.4 


Effect  of  Splitting  Error  on  32"''^  Fourier  Coefficient  with  Changing 
CFL  for  Supersonic  Flow 
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Figure  B.6  Effect  of  Splitting  Error  on  15*^  Fourier  Coefficient  with  Changing  CFL 
for  Subsonic  Flow 
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Figure  B.9  Effect  of  Splitting  Error  on  45*^  Fourier  Coefficient  with  Changing  CFL 
for  Subsonic  Flow 
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Appendix  C.  Determination  of  Optimal  Frequeney  Augmentation 

Inerements 

The  number  of  Fourier  frequencies  N  included  in  an  adapted  split-domain  harmonic 
balance  solution  can  have  a  large  influence  on  the  time  needed  to  compute  the  so¬ 
lution.  To  hnd  the  values  of  N  that  give  the  best  overall  performance,  a  numerical 
experiment  was  conducted.  In  this  experiment,  a  one-dimensional  split-domain  Eu¬ 
ler  solver  (see  Section  4.2)  was  run  a  fixed  number  of  iterations  for  every  number 
of  frequencies  N  from  3  to  100.  The  execution  time  for  each  run  was  recorded, 
normalized  by  the  execution  time  for  N  =  3,  and  plotted  in  Fig.  C.l. 


Figure  C.l  Effect  of  Series  Length  on  Split  Domain  Solver  Run  Time 

Values  of  N  that  result  in  locally  minimized  execution  time  were  identihed  as 
candidates  for  frequency  augmentation.  These  are  listed  in  Table  C.l.  Note  that 
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Table  C.l  Candidate  Numbers  of  Fourier  Frequencies,  N,  for  an  Efficient  Adaptive 
Harmonic  Balance  Solution 


0 

1 

2 
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4 

5 

6 
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10 

12 

13 

16 
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19 
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24 
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31 

32 

37 

38 

40 

45 

49 

52 

58 

62 

67 

73 

82 

87 

94 

97 

frequency  counts  0,  1,  and  2  are  also  included  in  this  table,  though  they  were  not 
included  in  the  experiment. 

The  only  computations  in  the  split-domain  harmonic  balance  method  that  are 
not  linearly  dependent  on  N  are  the  FFT  and  inverse  FFT  used  to  transform  the 
solution  vector  between  the  frequency  domain  and  the  time  domain.  The  values  of 
N  identihed  in  Table  C.l  are  therefore  those  that  minimized  the  cost  of  the  FFT 
calculations.  In  the  tested  solver  implementation,  FFT  calculations  were  performed 
via  “Fastest  Fourier  Transform  in  the  West”  (FFTW)  library  calls  (40).  The  locally 
minimizing  values  of  N  could  change  if  a  different  FFT  library  was  employed. 

While  any  adaptive  harmonic  balance  computation  could  include  every  value 
of  N  listed  in  Table  C.l,  some  values  could  be  omitted  if  information  is  known 
about  the  flow  being  calculated.  If  the  flow  held  is  known  to  contain  strong  moving 
discontinuities  throughout  the  entire  computational  domain,  then  it  will  have  a  high 
frequency  content  everywhere  as  well.  In  such  a  case,  it  is  unnecessary  to  include 
every  possible  low-frequency  candidate.  For  such  a  calculation,  the  following  values 
might  be  appropriate:  6,  12,  17,  24,  32,  37,  40,  45,  49,  ... .  In  this  set  of  N,  some 
lower  values  have  been  omitted,  because  it  is  known  that  the  lower  N  will  not  appear 
in  the  hnal  frequency  distribution.  The  low  N  that  are  included  simply  provide 
“stepping  stones”  to  the  hnal  solution.  All  of  the  higher  N  have  been  included, 
however,  so  that  the  hnal  frequency  distribution  can  be  ht  to  the  how  as  closely 
as  possible.  If,  on  the  other  hand,  the  how  held  being  modeled  contains  regions  of 
smoothly  unsteady  (or  even  steady-state)  how,  it  would  be  appropriate  to  include  all 
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of  the  candidate  N.  In  any  case,  it  will  not  affect  the  hnal  solution  if  all  candidate 
frequencies  are  included  in  a  calculation.  It  simply  may  take  longer  to  compute  than 
would  otherwise  be  necessary. 


C-3 


Appendix  D.  Implementation  of  Boundary  Conditions  for  the  1-D 

and  Quasi- 1-D  Euler  Equation 

Boundary  conditions  for  the  adaptive  split-domain  1-D  Euler  and  quasi- 1-D  Euler 
solvers  were  enforced  by  calculating  appropriate  values  for  a  single  ghost  cell  outside 
the  boundary  of  the  grid.  Ghost  cell  values  were  computed  after  the  first  integration 
of  the  frequency- domain  ODE  and  subsequent  transformation  to  the  time  domain, 
but  before  integration  of  the  time  domain  equation  (between  steps  2  and  3  of  the 
split-domain  update  process  as  dehned  in  Section  2.4,  page  2-10.)  This  approach 
allowed  the  use  of  conventional  steady-state  CFD  boundary  conditions. 

Supersonic  Boundary  Conditions 

A  characteristic  analysis  of  1-D  supersonic  flow  shows  that  all  flow  informa¬ 
tion  propagates  in  the  downstream  direction  (34).  Therefore,  at  a  supersonic  inflow 
boundary,  all  flow  information  comes  from  outside  the  computational  domain.  Con¬ 
versely,  at  a  supersonic  outflow  boundary  all  flow  information  comes  from  inside  the 
computational  domain. 

Supersonic  Inflow.  At  a  supersonic  inflow  boundary,  all  ghost  cell  values  must 
be  specihed.  At  the  inlet,  Mach  number,  density  and  pressure  were  given.  Inlet 
velocity  was  obtained  using  the  dehnition  of  Mach  number  and  the  relation  for  the 
speed  of  sound  in  a  perfect  gas: 


Unsteady  supersonic  boundary  conditions  were  obtained  by  varying  Mach  num¬ 
ber  while  holding  pressure  and  density  constant.  The  unsteady  Mach  number  was 
given  by 

M{t)  =  2.0  -|-  asin{ut)  (D.2) 
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where  a  and  u  were  the  amplitude  and  frequency  of  the  unsteady  variation.  For 
harmonic  balance  solutions,  samples  of  the  unsteady  Mach  number  were  obtained 
at  2N  +  1  points  spanning  one  period  of  oscillation,  where  N  was  the  number  of 
frequencies  included  in  the  harmonic  balance  solution  in  the  interior  cell  adjacent  to 
the  boundary. 


M{n)  =  2.0  +  a.  sin 


27rn  \ 
2N  +  l) 


n  =  0,l,...,2iV  (D.3) 


Conservative  flow  variables  were  constructed  from  the  specihed  primitive  vari¬ 


ables. 


Qg 


P 

pu 

f  (h%  +  b"). 


(D.4) 


Supersonic  Outflow.  At  a  supersonic  outflow  boundary,  ghost  cell  values  were 
obtained  by  means  of  a  zeroth-order  extrapolation  from  the  interior  cell  adjacent  to 
the  boundary. 
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Piston  Boundary  Condition 

The  action  of  an  oscillating  piston  was  simulated  by  specifying  a  sinusoidal 
velocity  in  the  ghost  point,  while  maintaining  a  zero  pressure  gradient  and  constant 
stagnation  enthalpy.  A  zero-pressure  gradient  was  maintained  by  setting  the  pressure 
in  the  ghost  cell  equal  to  the  pressure  in  the  interior  cell  adjacent  to  the  boundary, 
i.e.,  p  =  Padj-  For  inviscid  flow  of  a  perfect  gas,  stagnation  enthalpy,  given  by 
H  =  e  +  p/ p+^u^.,  is  constant  (41).  Therefore,  it  could  be  computed  from  the  initial 
conditions  in  the  tube.  Using  the  perfect  gas  relation  p  =  (7  —  l)pe,  and  the  fact 
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that  the  initial  velocity  was  zero,  the  stagnation  enthalpy  was  given  by 


H 


Pinit 

(7-  l)pMt' 


(D,6) 


For  time-accurate  solutions,  the  velocity  in  the  ghost  cell  was  given  by 


u{t)  =  asin(ci;t). 


(D,7) 


For  harmonic  balance  solutions,  samples  of  the  sinusoidal  velocity  were  obtained 
at  2N  +  1  points  spanning  one  period  of  oscillation,  where  N  was  the  number  of 
frequencies  included  in  the  harmonic  balance  solution  in  the  interior  cell  adjacent  to 
the  boundary. 


u{n)  =  a.  sin 


27rn  \ 
2N  +  1} 


n  =  0,l,...,2iV  (D.8) 


Given  m,  p,  and  H,  the  conservative  variables  in  the  ghost  cell  were  computed. 
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Subsonic  Characteristic  Farfield  Boundary  Condition 

A  characteristic  analysis  of  1-D  subsonic  flow  reveals  that  there  are  two  charac¬ 
teristics  traveling  downstream,  and  one  characteristic  traveling  upstream  (34).  Thus 
at  a  subsonic  inflow  boundary,  two  properties  must  be  specihed  and  one  must  be 
computed  from  the  interior  flow  held.  Conversely,  at  a  subsonic  outhow  bound¬ 
ary,  only  one  property  can  be  specihed,  and  two  properties  must  be  computed  from 
the  interior  how  held.  The  boundary  conditions  implemented  in  this  research  are 
developed  in  (34).  Only  the  implementation  is  presented  here. 
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Inflow.  Given  specified  farfield  properties  pf,  uj,  and  pf,  and  properties 
in  the  first  interior  cell,  padj,  u^dj,  and  padj,  define 


do  -  2^Pf  +  Ps^di) 

(D.IO) 

do  =  ^(d/+dadj) 

(D.ll) 

hpo 

Qq  —  .  - 

(D.12) 

V  do 

The  primitive  variables  at  the  boundary  are  given  by 

Pb  ~  Pf  T  dadj  T  PoO'oi'^f  ^adj) 

(D.13) 

,  Pf-Pb 

Pb  -  Pf  +  2 

(Iq 

(D.14) 

,  Pf-Pb 

Ub  =  Uf+— - 

Potto 

(D.15) 

The  primitive  variables  in  the  ghost  cells  were  obtained  through  a  hrst-order 
extrapolation  of  the  boundary  values  and  the  values  in  the  first  cell  in  the  grid. 


Pg 

‘^Pb  dadj 

(D.16) 

Ug 

=  ‘2Ub  Madj 

(D.17) 

Pg 

=  2pb  -  PadJ 

(D.18) 

These  primitive  variables  were  used  to  compute  conservative  variables  in  the  ghost 
point  using  Eq.  D.4. 

Outflow.  For  the  subsonic  outflow  boundary  condition,  pressure  was 
specihed,  and  density  and  velocity  were  computed.  The  primitive  variables  at  a 
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subsonic  outflow  boundary  were  given  by 


Pb 

=  Pf 

(D.19) 

pb 

Pb  Padj 

=  Padj  +  2 

Uq 

(D.20) 

Ub 

Padj  ~  Pb 

=  Madj  + 

PoQq 

(D.21) 

These  values  were  extrapolated  to  the  ghost  point  and  conservative  variables  were 
computed  using  Eq.  D.4. 

Subsonic  Nozzle  Outflow  Boundary  Condition 

A  variation  of  the  characteristic  subsonic  outflow  BC  was  incorporated  into 
the  quasi- 1-D  Euler  solver.  This  boundary  condition,  developed  in  (35),  takes  into 
account  the  varying  geometry  of  the  nozzle.  In  addition,  the  specihed  property  at 
the  boundary  is  density  rather  than  pressure. 

At  the  exit  boundary,  the  time-accurate  unsteady  density  was  specihed  accord¬ 
ing  to 

p{t)  =  Pavg  +  asin(a;t).  (D.22) 

For  harmonic  balance  solutions,  samples  of  the  unsteady  density  were  obtained  at 
2N  +  1  points  spanning  one  period  of  oscillation,  where  N  was  the  number  of  fre¬ 
quencies  included  in  the  harmonic  balance  solntion  in  the  interior  cell  adjacent  to 
the  bonndary. 

f  ^TTTl  \ 

p{n)  =  Pavg  +  asin  (  ^  n  =  0, 1, . . . ,  2iV  (D.23) 
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The  remaining  primitive  variables  at  the  exit  were  given  by  the  solution  to  two 


differential  equations  (35) 


(a  +  u) 


4  + 


f  dp  du\  1  dA  2 


0 

0 


(D.24a) 

(D.24b) 


In  Eq.  D.24b,  A  is  the  cross  sectional  area  of  the  nozzle  at  the  exit,  and  a  is  the  local 
speed  of  sound. 

Equation  D.24a  was  solved  for  Cb  by  approximating  the  spatial  derivatives  with 
hrst-order  backward  differences,  and  substituting  the  known  exit  density. 


eb 


Cadj 

l-(7-l)(l-Padj/p) 
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The  exit  velocity,  u,  was  obtained  by  solving  Eq.  D.24b.  The  spatial  derivatives 
were  replaced  with  hrst-order  backward  differences  and  terms  were  rearranged  to 
form  a  quadratic  expression  for  u.  Of  the  two  possible  solutions,  one  produced  a 
non-physical  answer  when  the  nozzle  cross  sectional  area  was  assumed  constant  and 
the  pressure  gradient  ||  was  negative.  This  solution  was  discarded,  leaving  the 
following  equation  for  u: 


-P  +  \/P^  +  4pa^(padj  -  P  +  pau^dj) 
2  pa 


(D.26) 


where 


p  =  pa^- - h  (p  -  Padj)  +  pa{a  -  Uadj)  (D.27) 

Tiexit 

and  the  quantity  AA  was  the  change  in  nozzle  cross  sectional  area  across  the  cell 
adjacent  to  the  boundary. 
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The  quantities  r,  u,  and  e,  were  used  to  compute  conservative  variables  in  the 
ghost  point  using  Eq.  D.4. 
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